Miscellaneous functions¶

Return the adjugate and the determinant of a square integer matrix.

The adjugate is the transpose of the cofactor matrix. It satisfies \(A \operatorname{adj}(A) = \operatorname{adj}(A) A = \det(A) I\), and it is defined even when \(A\) is singular. All computations are done exactly with Python integers, so the result is correct at any magnitude.

Parameters:

• mat (ArrayLike) –

  A square integer matrix.

Returns:

• tuple[ndarray, int] –

  The adjugate and the determinant. The adjugate has int64 dtype when its entries fit, and object dtype (Python integers) otherwise.

Raises:

• ValueError –

  If the matrix is not square.

Examples:

We compute the adjugate and the determinant of a \(2\times 2\) matrix.

>>> from cytools.utils import adjugate
>>> adjugate([[1, 2], [3, 4]])
(array([[ 4, -2],
       [-3,  1]]), -2)

Source on GitHub

Convert a float array to an array of fmpq, FLINT's rational type.

This is a shortcut for array_to_flint with t=float. Each entry is converted with float_to_fmpq.

Parameters:

• arr (ArrayLike) –

  An array with float entries.

Returns:

• ndarray –

  An array of the same shape, with object dtype and fmpq entries.

Examples:

We convert a float array to an fmpq array.

>>> from cytools.utils import array_float_to_fmpq
>>> array_float_to_fmpq([0.5, 0.25])
array([1/2, 1/4], dtype=object)

Source on GitHub

Convert an array of fmpq, FLINT's rational type, to a float array.

This is a shortcut for array_from_flint with t=flint.fmpq.

Parameters:

• arr (ndarray) –

  An array with fmpq entries.

Returns:

• ndarray –

  An array of the same shape, with float dtype.

Source on GitHub

Convert an array of fmpz, FLINT's integer type, to an integer array.

This is a shortcut for array_from_flint with t=flint.fmpz.

Parameters:

• arr (ndarray) –

  An array with fmpz entries.

Returns:

• ndarray –

  An array of the same shape, with int dtype.

Source on GitHub

Convert an array of FLINT numbers to an array of integers or floats.

Entries of type fmpz (FLINT's integer type) are converted to integers, and entries of type fmpq (FLINT's rational type) to floats.

Parameters:

• arr (ndarray) –

  An array with fmpz or fmpq entries.

• t (type | None, default: None ) –

  The type of the entries, either flint.fmpz or flint.fmpq. If not specified, it is taken to be the type of the first entry of arr.

Returns:

• ndarray –

  An array of the same shape, with int or float dtype.

Raises:

• ValueError –

  If the entries are not of type fmpz or fmpq.

Examples:

We convert an fmpq array to a float array.

>>> import numpy as np
>>> from flint import fmpq
>>> from cytools.utils import array_from_flint
>>> array_from_flint(np.array([fmpq(1, 2), fmpq(1, 3)]))
array([0.5       , 0.33333333])

Source on GitHub

Convert an integer array to an array of fmpz, FLINT's integer type.

This is a shortcut for array_to_flint with t=int.

Parameters:

• arr (ArrayLike) –

  An array with integer entries.

Returns:

• ndarray –

  An array of the same shape, with object dtype and fmpz entries.

Examples:

We convert an integer array to an fmpz array.

>>> from cytools.utils import array_int_to_fmpz
>>> array_int_to_fmpz([[1, 0], [0, 2]])
array([[1, 0],
       [0, 2]], dtype=object)

Source on GitHub

Convert a numpy array to an array of FLINT numbers.

Depending on t, each entry is converted either to an fmpz (FLINT's integer type) or to an fmpq (FLINT's rational type), the latter using float_to_fmpq. The functions array_int_to_fmpz and array_float_to_fmpq are shortcuts for the two cases.

Parameters:

• arr (ndarray) –

  The array to convert.

• t (type | None, default: None ) –

  The type of the entries of arr. If it is int, the entries are converted to fmpz, and otherwise to fmpq. If not specified, arr.dtype is used. Since a numpy dtype is never int, this always results in a conversion to fmpq.

Returns:

• ndarray –

  An array of the same shape, with object dtype and fmpz or fmpq entries.

Examples:

We convert an integer array and a float array.

>>> import numpy as np
>>> from cytools.utils import array_to_flint
>>> array_to_flint(np.array([[1, 0], [0, 2]]), t=int)
array([[1, 0],
       [0, 2]], dtype=object)
>>> array_to_flint(np.array([0.5, 0.25]))
array([1/2, 1/4], dtype=object)

Source on GitHub

Fetch reflexive polytopes from the Kreuzer-Skarke or Schöller-Skarke database.

Four-dimensional polytopes are fetched from the Kreuzer-Skarke database, and five-dimensional ones from the Schöller-Skarke database. The polytopes can be filtered by the Hodge numbers and Euler characteristic of the associated Calabi-Yau hypersurfaces, and, for four-dimensional polytopes, by their numbers of points, vertices, and facets.

This function can also be imported from the package root.

Parameters:

• h11 (int | None, default: None ) –

  The Hodge number \(h^{1,1}\) of the Calabi-Yau hypersurface.

• h12 (int | None, default: None ) –

  The Hodge number \(h^{1,2}\) of the Calabi-Yau hypersurface.

• h13 (int | None, default: None ) –

  The Hodge number \(h^{1,3}\) of the Calabi-Yau hypersurface. Only used for five-dimensional polytopes.

• h21 (int | None, default: None ) –

  The Hodge number \(h^{2,1}\) of the Calabi-Yau hypersurface. This is equivalent to h12.

• h22 (int | None, default: None ) –

  The Hodge number \(h^{2,2}\) of the Calabi-Yau hypersurface. Only used for five-dimensional polytopes.

• h31 (int | None, default: None ) –

  The Hodge number \(h^{3,1}\) of the Calabi-Yau hypersurface. This is equivalent to h13.

• chi (int | None, default: None ) –

  The Euler characteristic of the Calabi-Yau hypersurface.

• lattice (str | None, default: 'N' ) –

  The lattice in which the polytopes are defined, "N" or "M". The Hodge numbers and the Euler characteristic refer to the Calabi-Yau hypersurface associated with polytopes in this lattice, and favorability is checked with respect to it. It must be specified if any Hodge number, the Euler characteristic, or favorable is specified.

• dim (int, default: 4 ) –

  The dimension of the polytopes. The options are 4 and 5.

• n_points (int | None, default: None ) –

  The number of lattice points of the polytopes. Only used for four-dimensional polytopes.

• n_vertices (int | None, default: None ) –

  The number of vertices of the polytopes. Only used for four-dimensional polytopes.

• n_dual_points (int | None, default: None ) –

  The number of lattice points of the dual polytopes. Only used for four-dimensional polytopes.

• n_facets (int | None, default: None ) –

  The number of facets of the polytopes. Only used for four-dimensional polytopes.

• limit (int, default: 1000 ) –

  The maximum number of polytopes to return.

• samples (int | None, default: None ) –

  If specified, this many polytopes are chosen at random, without replacement, among the (at most limit) fetched ones. It must not be larger than limit.

• sample_seed (int | None, default: None ) –

  The seed of the random number generator used for sampling. If not specified, a different sample is drawn every time.

• timeout (int, default: 60 ) –

  The maximum number of seconds to wait for the server to return the data.

• as_list (bool, default: True ) –

  Whether to return a list of polytopes instead of a generator.

• backend (str | None, default: None ) –

  The backend used to construct the convex hull of each Polytope. The options are "ppl", "qhull", and "palp". If not specified, it is chosen automatically.

• deterministic_glsm_basis (bool, default: False ) –

  Whether to choose the GLSM (and thus divisor) basis of each polytope in a deterministic manner, so that it is consistent across different machines. The basis chosen with this flag may differ from the default one.

• dualize (bool, default: False ) –

  Whether to return the dual of each polytope instead of the polytope itself.

• favorable (bool | None, default: None ) –

  If true, only favorable polytopes are returned, and if false, only non-favorable ones. Favorability is checked before dualizing. If not specified, all polytopes are returned.

• verbosity (int, default: 0 ) –

  The verbosity level. If it is 1 or larger, the URL of the request is printed, and if it is 2 or larger, the received data is printed as well.

Returns:

• Generator[Polytope, None, None] | list[Polytope] –

  A list of Polytope objects, or a generator of them if as_list is false.

Raises:

• ValueError –

  If dim or lattice is not one of the allowed options, if lattice is not specified when required, if samples is larger than limit, if both h12 and h21 (or h13 and h31) are given with different values, or if the given Hodge numbers and Euler characteristic are inconsistent.

Examples:

We fetch polytopes from the Kreuzer-Skarke and Schöller-Skarke databases with a few different parameters.

>>> from cytools import fetch_polytopes
>>> g = fetch_polytopes(h11=27, as_list=False)  # a generator of polytopes
>>> print(next(g))
A 4-dimensional reflexive lattice polytope in ZZ^4
>>> polys = fetch_polytopes(h11=27, limit=100)  # a list of polytopes
>>> len(polys)
100
>>> g_5d = fetch_polytopes(h11=1000, dim=5, limit=100, as_list=False)
>>> print(next(g_5d))
A 5-dimensional reflexive lattice polytope in ZZ^5

Source on GitHub

Restrict a sparse tensor to a subset of its indices.

Only the entries whose indices all belong to indices are kept, and they are reindexed so that indices[i] becomes i. The indices of each resulting key are sorted in ascending order. This can be used, for example, to obtain the intersection numbers in a basis of divisors that is a subset of the prime toric divisors.

Parameters:

• tensor (dict) –

  A symmetric sparse tensor, as a dictionary of the form {(a, b, ..., c): M_ab...c, ...}.

• indices (list[int]) –

  The indices to keep. Their order determines the new indexing.

Returns:

• dict –

  The filtered and reindexed sparse tensor, in the same format as the input.

Examples:

We construct a simple tensor and keep only the indices 1 and 3.

>>> from cytools.utils import filter_tensor_indices
>>> tensor = {(0,1):0, (1,1):1, (1,2):2, (1,3):3, (2,3):4}
>>> filter_tensor_indices(tensor, [1,3])
{(0, 0): 1, (0, 1): 3}

The intersection numbers of a toric variety in its divisor basis can be obtained in this way from those of all the prime toric divisors.

>>> from cytools import Polytope
>>> p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
>>> v = p.triangulate().get_toric_variety()
>>> intnums_all = v.intersection_numbers(in_basis=False)
>>> intnums_basis = v.intersection_numbers(in_basis=True)
>>> intnums_basis == filter_tensor_indices(intnums_all, v.divisor_basis())
True

Source on GitHub

Find new points that make a set of points full-dimensional.

If the input points affinely span a \(d\)-dimensional subspace of \(\mathbb{R}^n\), this function returns \(n-d\) integer points such that the input points together with the new ones affinely span all of \(\mathbb{R}^n\). This is useful to turn a point configuration that is not full-dimensional into one that is.

Parameters:

• pts (ArrayLike) –

  The points, one per row. There must be at least one.

Returns:

• ndarray –

  The \(n-d\) new points, one per row.

Raises:

• ValueError –

  If the list of points is empty.

Examples:

We construct three points spanning a plane in 3 dimensions, and find a new point off that plane.

>>> from cytools.utils import find_new_affinely_independent_points
>>> find_new_affinely_independent_points([[1,0,1],[0,0,1],[0,1,1]])
array([[1, 0, 2]])

Source on GitHub

Convert a float to an fmpq, FLINT's rational number type.

The float is approximated by the closest rational number with denominator at most \(10^6\) (see fractions.Fraction.limit_denominator), which recovers simple fractions despite floating-point rounding errors.

Parameters:

• c (float) –

  The number to convert.

Returns:

• fmpq –

  The rational approximation of c.

Examples:

We convert a few floats to rational numbers.

>>> from cytools.utils import float_to_fmpq
>>> float_to_fmpq(0.1), float_to_fmpq(0.333333333333), float_to_fmpq(2.45)
(1/10, 1/3, 49/20)

Source on GitHub

Convert an fmpq, FLINT's rational number type, to a float.

Parameters:

• c (fmpq) –

  The rational number to convert.

Returns:

• float –

  The number as a float.

Examples:

We convert a few rational numbers to floats.

>>> from flint import fmpq
>>> from cytools.utils import fmpq_to_float
>>> [fmpq_to_float(fmpq(p, q)) for p, q in [(1, 2), (1, 3), (49, 20)]]
[0.5, 0.3333333333333333, 2.45]

Source on GitHub

Return the greatest common divisor of two floats.

The gcd is computed with the Euclidean algorithm, where a remainder smaller than tol in absolute value is treated as zero. The result is thus, up to rounding errors, the largest number that divides both a and b a whole number of times. This is useful since the standard gcd functions do not accept non-integer values.

Parameters:

• a (float) –

  The first number.

• b (float) –

  The second number.

• tol (float, default: 1e-05 ) –

  The tolerance below which a remainder is considered to be zero.

Returns:

• float –

  The (non-negative) gcd of a and b.

Examples:

We compute the gcd of two floats. The exact answer is 0.1, but there are small rounding errors.

>>> from cytools.utils import gcd_float
>>> gcd_float(0.2, 0.5)
0.09999999999999998

Source on GitHub

Return the greatest common divisor of a list of numbers.

The gcd is computed pairwise with gcd_float, using its default tolerance, so the entries do not need to be integers.

Parameters:

• arr (Iterable[float]) –

  The numbers. It must not be empty.

Returns:

• float –

  The gcd of all the numbers.

Examples:

We compute the gcd of a few integers.

>>> from cytools.utils import gcd_list
>>> gcd_list([4, 6, 10])
2

Source on GitHub

Convert a height vector to the corresponding divisor class in the basis.

This is the same as project_heights_to_kahler, but only the \(h^{1,1}\) coefficients of the divisor class in the basis of divisors are returned.

This is one of a few temporary helpers for converting between heights and Kähler parameters. It assumes that poly is a four-dimensional reflexive polytope whose prime toric divisors are the \(h^{1,1}+4\) nonzero lattice points not interior to facets (i.e. that poly is favorable), and it uses the basis of divisors given by Polytope.glsm_basis.

Parameters:

• poly (Polytope) –

  The polytope.

• heights_in (ArrayLike) –

  The heights, a vector of length \(h^{1,1}+5\) with the height of the origin first, followed by those of the prime toric divisors.

• prime_divisors (ArrayLike | None, default: None ) –

  The prime toric divisors that are not in the basis, written as linear combinations of the basis divisors. See project_heights_to_kahler for details.

Returns:

• ndarray –

  The \(h^{1,1}\) coefficients of the divisor class \(\sum_i (h_i - h_0) D_i\) in the basis of divisors.

Source on GitHub

Return a decorator that caches the results of a method per instance.

This works like functools.lru_cache, but each instance gets its own cache, stored in its _cache dictionary (which is created if needed), so that cached results are never shared between instances. Resetting _cache to an empty dictionary clears the cache of the instance. The arguments of the decorated method must be hashable.

Parameters:

• maxsize (int | None, default: 128 ) –

  The maximum number of results cached per instance and method, as in functools.lru_cache. If None, the cache can grow without bound.

Returns:

• Callable –

  The decorator.

Examples:

We define a class with a cached method. The result is only computed the first time it is requested.

>>> from cytools.utils import instanced_lru_cache
>>> class A:
...     @instanced_lru_cache(maxsize=None)
...     def square(self, x):
...         print("computing")
...         return x**2
>>> a = A()
>>> a.square(3)
computing
9
>>> a.square(3)
9

Source on GitHub

Return a basis of the integral nullspace of an integer matrix.

The returned vectors form a basis of the lattice \(\{x \in \mathbb{Z}^n : Mx = 0\}\), where \(n\) is the number of columns of M. This lattice is saturated: if \(kx\) lies in it for some nonzero integer \(k\), then so does \(x\). As a consequence, the basis can always be completed to a basis of \(\mathbb{Z}^n\). Note that dividing the vectors of a rational nullspace basis by their gcds is not enough to obtain such a basis; in general they only span a sublattice of finite index. For example, for M = [[9, 4, 7]] the kernel vector [1, -4, 1] might not be an integer combination of them.

The basis is computed from the Hermite normal form \(H = UM^T\), where \(U\) is unimodular. The rows of \(U\) corresponding to zero rows of \(H\) span the kernel, and since \(U\) is unimodular they extend to a basis of \(\mathbb{Z}^n\).

Parameters:

• M (ArrayLike) –

  The integer matrix whose nullspace is computed.

• reduce_by_gcd (bool, default: True ) –

  Ignored. It is only accepted for backwards compatibility, since the returned basis vectors are always primitive.

Returns:

• ndarray –

  The basis vectors, as the columns of an array of shape (n, k), where k is the dimension of the nullspace. The entries have int64 dtype when they fit, and are Python integers (object dtype) otherwise.

Examples:

We compute the integral nullspace of a \(1\times 3\) matrix. The kernel vector [1, -4, 1] is 5 times the first column plus the second one.

>>> from cytools.utils import integral_nullspace
>>> integral_nullspace([[9, 4, 7]])
array([[ -4,  21],
       [  9, -49],
       [  0,   1]])

Source on GitHub

Convert a divisor class in the basis to a height vector.

The heights of the basis divisors are set to the given coefficients, and all the other heights, including that of the origin, are set to zero. This is a right inverse of heights_to_kahler.

This is one of a few temporary helpers for converting between heights and Kähler parameters. It assumes that poly is a four-dimensional reflexive polytope whose prime toric divisors are the \(h^{1,1}+4\) nonzero lattice points not interior to facets (i.e. that poly is favorable), and it uses the basis of divisors given by Polytope.glsm_basis.

Parameters:

• poly (Polytope) –

  The polytope.

• kahler_in (ArrayLike) –

  The \(h^{1,1}\) coefficients of the divisor class in the basis of divisors.

Returns:

• ndarray –

  The heights, a vector of length \(h^{1,1}+5\) with the height of the origin first, followed by those of the prime toric divisors.

Source on GitHub

Return the index of the lattice generated by the rows of an integer matrix.

The index is computed as the product of the invariant factors (the diagonal entries of the Smith normal form) of mat. For a square matrix of full rank, this is \(|\det(\text{mat})|\), the index of the generated sublattice in \(\mathbb{Z}^n\). More generally, if the rank of mat equals its number of rows or of columns, this is the index of the generated lattice in its saturation (the set of integer points in its real span). Otherwise, the result is 0.

Parameters:

• mat (ArrayLike) –

  An integer matrix whose rows generate the lattice.

Returns:

• int –

  The lattice index.

Examples:

The vectors \((1,1)\) and \((1,-1)\) generate the sublattice of index 2 of \(\mathbb{Z}^2\) consisting of the points with even coordinate sum.

>>> from cytools.utils import lattice_index
>>> lattice_index([[1, 1], [1, -1]])
2

Source on GitHub

Apply LLL reduction to a set of lattice points.

The points are the rows of pts_in. The LLL algorithm is applied to the rows of pts_in.T, i.e. the vectors of coordinates of the points. This amounts to a unimodular change of lattice coordinates \(p \mapsto Ap\), with \(A \in GL(n, \mathbb{Z})\), that makes these vectors short and nearly orthogonal. If the points span a sublattice of dimension \(d < n\), the coordinates that vanish for all points come first, so that the last \(d\) columns of the result give the coordinates of the points in the sublattice.

Parameters:

• pts_in (ArrayLike) –

  The integer points, one per row.

• transform (bool, default: False ) –

  Whether to also return the transformation matrix and its inverse.

Returns:

• ndarray | tuple[ndarray, tuple[ndarray, ndarray]] –

  The reduced points, one per row. If transform is true, a tuple (pts_red, (A, Ainv)) is returned instead, where A is the unimodular matrix satisfying pts_red.T == A @ pts_in.T and Ainv is its inverse.

Examples:

We reduce a set of points spanning a 2-dimensional sublattice of \(\mathbb{Z}^3\). The first coordinate of all the reduced points vanishes.

>>> from cytools.utils import lll_reduce
>>> lll_reduce([[1,1,0],[0,1,1],[1,2,1],[0,0,0]])
array([[0, 1, 0],
       [0, 0, 1],
       [0, 1, 1],
       [0, 0, 0]])

Source on GitHub

Read polytopes from a file or a string, yielding them one at a time.

The polytopes can be specified by their vertices, in the format used in the Kreuzer-Skarke database, or by weight systems.

In the "ks" format, each polytope consists of a header line containing M:, whose first two numbers are the dimensions n m of the matrix that follows, and n lines with m integers each. The vertices are the columns of this matrix if it has more columns than rows, and its rows otherwise. In the "ws" format, each non-empty line is a weight system, such as 5 1 1 1 1 1 (the degree followed by the weights), and the polytope is the one that PALP associates with it.

This function is not meant to be called directly. Instead, it is used by read_polytopes and fetch_polytopes.

Parameters:

• input (str) –

  The name of the file to read, or the string containing the polytopes, depending on input_type.

• input_type (str, default: 'file' ) –

  Whether input is the name of a file or a string containing the data. The options are "file" and "str".

• format (str, default: 'ks' ) –

  The format of the data. The options are "ks", the format used in the Kreuzer-Skarke database, where the polytopes are specified by their vertices, and "ws", where each line is a weight system.

• backend (str | None, default: None ) –

  The backend used to construct the convex hull of each Polytope. The options are "ppl", "qhull", and "palp". If not specified, it is chosen automatically.

• deterministic_glsm_basis (bool, default: False ) –

  Whether to choose the GLSM (and thus divisor) basis of each polytope in a deterministic manner, so that it is consistent across different machines. The basis chosen with this flag may differ from the default one.

• dualize (bool, default: False ) –

  Whether to return the dual of each polytope instead of the polytope itself.

• favorable (bool | None, default: None ) –

  If true, only favorable polytopes are returned, and if false, only non-favorable ones. Favorability is checked before dualizing. If not specified, all polytopes are returned.

• lattice (str | None, default: None ) –

  The lattice, "N" or "M", with respect to which favorability is checked. It must be specified if favorable is.

• limit (int | None, default: None ) –

  The maximum number of polytopes to return. If not specified, all of them are returned.

Yields:

• Polytope –

  The polytopes that were read.

Raises:

• ValueError –

  If favorable is specified without lattice, if input_type or format is not one of the allowed options, or if the dimensions of a vertex matrix do not match its header. Since this is a generator, these errors are raised when it is first advanced.

Examples:

We read a polytope from a string taken from the Kreuzer-Skarke database.

>>> from cytools.utils import polytope_generator
>>> poly_data = '''4 5  M:10 5 N:376 5 H:272,2 [540]
...    1    0    0    0   -9
...    0    1    0    0   -6
...    0    0    1    0   -1
...    0    0    0    1   -1
... '''
>>> g = polytope_generator(poly_data, input_type="str")
>>> print(next(g))
A 4-dimensional reflexive lattice polytope in ZZ^4

Source on GitHub

Project a height vector onto the basis of divisors.

Given heights \(h_0, h_1, \ldots\) for the origin and the prime toric divisors \(D_1, D_2, \ldots\) of the polytope, this computes the divisor class \(\sum_i (h_i - h_0) D_i\), written in the basis of divisors. The prime toric divisors that are not in the basis are eliminated using prime_divisors, which gives each of them as a linear combination of the basis divisors.

This is one of a few temporary helpers for converting between heights and Kähler parameters. It assumes that poly is a four-dimensional reflexive polytope whose prime toric divisors are the \(h^{1,1}+4\) nonzero lattice points not interior to facets (i.e. that poly is favorable), and it uses the basis of divisors given by Polytope.glsm_basis.

Parameters:

• poly (Polytope) –

  The polytope.

• heights_in (ArrayLike) –

  The heights, a vector of length \(h^{1,1}+5\) with the height of the origin first, followed by those of the prime toric divisors.

• prime_divisors (ArrayLike | None, default: None ) –

  The prime toric divisors that are not in the basis, in order, written as linear combinations of the basis divisors, one per row. If not specified, they are read off from the rays of the toric effective cone of a Calabi-Yau hypersurface obtained from a triangulation of poly.

Returns:

• ndarray –

  A vector of length \(h^{1,1}+5\), indexed like heights_in, whose entries at the positions of the basis divisors are the coefficients of the divisor class, and whose other entries are zero.

Source on GitHub

Read polytopes from a file or a string.

The polytopes can be specified by their vertices, in the format used in the Kreuzer-Skarke database, or by weight systems.

In the "ks" format, each polytope consists of a header line containing M:, whose first two numbers are the dimensions n m of the matrix that follows, and n lines with m integers each. The vertices are the columns of this matrix if it has more columns than rows, and its rows otherwise. In the "ws" format, each non-empty line is a weight system, such as 5 1 1 1 1 1 (the degree followed by the weights), and the polytope is the one that PALP associates with it.

This function can also be imported from the package root.

Parameters:

• input (str) –

  The name of the file to read, or the string containing the polytopes, depending on input_type.

• input_type (str, default: 'file' ) –

  Whether input is the name of a file or a string containing the data. The options are "file" and "str".

• format (str, default: 'ks' ) –

  The format of the data. The options are "ks", the format used in the Kreuzer-Skarke database, where the polytopes are specified by their vertices, and "ws", where each line is a weight system.

• backend (str | None, default: None ) –

  The backend used to construct the convex hull of each Polytope. The options are "ppl", "qhull", and "palp". If not specified, it is chosen automatically.

• deterministic_glsm_basis (bool, default: False ) –

  Whether to choose the GLSM (and thus divisor) basis of each polytope in a deterministic manner, so that it is consistent across different machines. The basis chosen with this flag may differ from the default one.

• as_list (bool, default: False ) –

  Whether to return a list of polytopes instead of a generator.

• dualize (bool, default: False ) –

  Whether to return the dual of each polytope instead of the polytope itself.

• favorable (bool | None, default: None ) –

  If true, only favorable polytopes are returned, and if false, only non-favorable ones. Favorability is checked before dualizing. If not specified, all polytopes are returned.

• lattice (str | None, default: None ) –

  The lattice, "N" or "M", with respect to which favorability is checked. It must be specified if favorable is.

• limit (int | None, default: None ) –

  The maximum number of polytopes to return. If not specified, all of them are returned.

Returns:

• Generator[Polytope, None, None] | list[Polytope] –

  A generator of Polytope objects, or a list of them if as_list is true.

Raises:

• ValueError –

  If favorable is specified without lattice, if input_type or format is not one of the allowed options, or if the dimensions of a vertex matrix do not match its header. If a generator is returned, these errors are raised when it is first advanced.

Examples:

We read a polytope from a string taken from the Kreuzer-Skarke database.

>>> from cytools import read_polytopes
>>> poly_data = '''4 5  M:10 5 N:376 5 H:272,2 [540]
...    1    0    0    0   -9
...    0    1    0    0   -6
...    0    0    1    0   -1
...    0    0    0    1   -1
... '''
>>> polys = read_polytopes(poly_data, input_type="str", as_list=True)
>>> len(polys)
1
>>> print(polys[0])
A 4-dimensional reflexive lattice polytope in ZZ^4
>>> polys[0].h11(lattice="N"), polys[0].h21(lattice="N")
(2, 272)

We can also construct polytopes from weight systems. Here we use the weight system of the quintic, and take the dual of the resulting polytope.

>>> ws_polys = read_polytopes("5 1 1 1 1 1", input_type="str", format="ws",
...                           dualize=True, as_list=True)
>>> ws_polys[0].h11(lattice="N")
1

Source on GitHub

Set the basis of curves of a toric variety or Calabi-Yau.

Specifying a basis of curves also determines the dual basis of divisors. The basis can be given either as a vector of indices, or as a matrix whose rows are the basis curves, specified by their intersection numbers with the prime toric divisors. When a vector is used, this is equivalent to calling set_divisor_basis with the same vector. All cached quantities that depend on the basis are cleared.

This function is not meant to be called directly. Instead, it is called by the methods set_curve_basis of ToricVariety and set_curve_basis of CalabiYau. Generic bases given by matrices are an experimental feature, which must first be enabled with enable_experimental_features.

Parameters:

• tv_or_cy (ToricVariety | CalabiYau) –

  The toric variety or Calabi-Yau whose basis is set.

• basis (ArrayLike) –

  The basis. If it is a vector, its entries are the indices of the elements of the standard basis of the dual to the lattice of prime toric divisors, i.e. the curve basis dual to the divisor basis with the same indices. If it is a matrix, each row is a basis curve, given by its intersection numbers with the prime toric divisors, with or without a first column for the origin (i.e. the canonical divisor), as inferred from its shape. If the first column is omitted, it is set to minus the sum of the other entries of the row.

• include_origin (bool, default: True ) –

  Whether the indices of a vector basis count the origin as index 0. If false, the indices are shifted by one. Ignored for matrix bases.

Raises:

• ValueError –

  If the vector basis is invalid (see set_divisor_basis), if the matrix has the wrong rank or shape, if its rows are not valid curve charge vectors (i.e. the corresponding combinations of the points of the polytope, each with a 1 appended, do not vanish), or if they do not form an integral basis.

• Exception –

  If a matrix basis is given without enabling the experimental features.

Examples:

This function is not meant to be used directly, but it is called in the following example. We consider a toric variety with two independent curves, find the default basis of curves, and then set a basis of our choice.

>>> from cytools import Polytope
>>> p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
>>> t = p.triangulate()
>>> v = t.get_toric_variety()
>>> v.curve_basis()  # the default basis
array([5, 6])
>>> v.set_curve_basis([1,6])  # here we set a basis
>>> v.curve_basis()  # we get the basis we set
array([1, 6])
>>> v.curve_basis(as_matrix=True)  # the basis in matrix form
array([[-18,   1,   9,   6,   1,   1,   0],
       [ -6,   0,   3,   2,   0,   0,   1]])

Since a vector is used, this is the same as calling set_divisor_basis with the same vector. The two functions only differ for generic bases given by matrices.

Source on GitHub

Set the basis of divisors of a toric variety or Calabi-Yau.

Specifying a basis of divisors also determines the dual basis of curves. The basis can be given either as a vector of indices of prime toric divisors, or as a matrix whose rows are linear combinations of prime toric divisors. All cached quantities that depend on the basis are cleared.

This function is not meant to be called directly. Instead, it is called by the methods set_divisor_basis of ToricVariety and set_divisor_basis of CalabiYau. Generic bases given by matrices are an experimental feature, which must first be enabled with enable_experimental_features.

Parameters:

• tv_or_cy (ToricVariety | CalabiYau) –

  The toric variety or Calabi-Yau whose basis is set.

• basis (ArrayLike) –

  The basis. If it is a vector, its entries are the indices of the prime toric divisors (i.e. of the corresponding points of the polytope) that form the basis. The indices are sorted, so their order does not matter. If it is a matrix, each row is a basis divisor written as a linear combination of the prime toric divisors, with or without a first column for the origin (i.e. the canonical divisor), as inferred from its shape.

• include_origin (bool, default: True ) –

  Whether the indices of a vector basis count the origin as index 0. If false, the indices are shifted by one. Ignored for matrix bases.

Raises:

• ValueError –

  If the indices are out of range, if the matrix has the wrong rank or shape, or if the given divisors do not form an integral basis.

• Exception –

  If a matrix basis is given without enabling the experimental features.

Examples:

This function is not meant to be used directly, but it is called in the following example. We consider a toric variety with two independent divisors, find the default basis, and then set a basis of our choice.

>>> from cytools import Polytope
>>> p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
>>> t = p.triangulate()
>>> v = t.get_toric_variety()
>>> v.divisor_basis()  # the default basis
array([5, 6])
>>> v.set_divisor_basis([1,6])  # here we set a basis
>>> v.divisor_basis()  # we get the basis we set
array([1, 6])
>>> v.divisor_basis(as_matrix=True)  # the basis in matrix form
array([[0, 1, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 1]])

Source on GitHub

Solve the sparse linear system \(Mx + C = 0\).

The system is solved in the least-squares sense, through the normal equations \(M^T M x = -M^T C\). The "sksparse" backend uses a sparse Cholesky factorization from scikit-sparse (an optional dependency), and the "scipy" backend uses scipy.sparse.linalg.spsolve. If check is true, a solution is only accepted if every entry of the residual \(Mx + C\) is at most backend_error_tol in absolute value.

Parameters:

• M (csr_matrix) –

  The coefficient matrix, as a SciPy sparse matrix.

• C (list[float]) –

  The constant term.

• backend (str, default: 'all' ) –

  The solver to use. The options are "all", "sksparse", and "scipy". With "all", the "sksparse" and "scipy" backends are tried in that order until one of them finds a solution.

• check (bool, default: True ) –

  Whether to check the solution by computing the residual.

• backend_error_tol (float, default: 0.0001 ) –

  The largest absolute value of the residual entries for which the solution is accepted. Only used if check is true.

• verbosity (int, default: 0 ) –

  The verbosity level. If it is 1 or larger, a message is printed whenever a backend fails.

Returns:

• list[float] | ndarray | None –

  The solution, or None if no acceptable solution was found. The "scipy" backend returns a list, and the "sksparse" backend returns an array.

Raises:

• ValueError –

  If backend is not one of the allowed options.

• RuntimeError –

  If backend is "all" and every backend raised an error, rather than just failing to find a solution.

Examples:

We solve a very simple linear system.

>>> from cytools.utils import to_sparse, solve_linear_system
>>> id_array = [[0,0,1],[1,1,1],[2,2,1],[3,3,1],[4,4,1]]
>>> M = to_sparse(id_array, sparse_type="csr")
>>> C = [1,1,1,1,1]
>>> solve_linear_system(M, C)
[-1.0, -1.0, -1.0, -1.0, -1.0]

Source on GitHub

Convert a dense symmetric tensor to a sparse tensor.

The output is a dictionary {(a, b, ..., c): M_ab...c, ...} containing only the nonzero entries with sorted indices \(a \le b \le \cdots \le c\). Optionally, a change of basis is applied first, in the same way as in symmetric_sparse_to_dense.

Parameters:

• tensor (ArrayLike) –

  A dense symmetric tensor, with the same length along all of its axes.

• basis (ArrayLike | None, default: None ) –

  A matrix whose rows are the new basis elements, written in terms of the original ones. Currently, only square matrices are handled correctly. If not specified, no change of basis is applied.

Returns:

• dict –

  The sparse tensor, as a dictionary of the form {(a, b, ..., c): M_ab...c, ...}. The values are numpy scalars.

Raises:

• ValueError –

  If the axes of tensor do not all have the same length.

Examples:

We convert a dense symmetric matrix to a sparse one.

>>> from cytools.utils import symmetric_dense_to_sparse
>>> symmetric_dense_to_sparse([[1,2],[2,3]])
{(0, 0): np.int64(1), (0, 1): np.int64(2), (1, 1): np.int64(3)}

Source on GitHub

Convert a symmetric sparse tensor to a dense array.

The sparse tensor is a dictionary {(a, b, ..., c): M_ab...c, ...}, where each entry only needs to be given for one ordering of its indices; the dense tensor is filled in for all their permutations. Optionally, a change of basis is then applied, so that the entries of the output are $$ T'{ij\cdots k} = \sum\, T_{ab\cdots c}\,, $$ where } B_{ia} B_{jb} \cdots B_{kc\(B\) is the basis matrix.

Parameters:

• tensor (dict) –

  The symmetric sparse tensor, as a dictionary of the form {(a, b, ..., c): M_ab...c, ...}.

• basis (ArrayLike | None, default: None ) –

  A matrix whose rows are the new basis elements, written in terms of the original ones. It must have at least as many columns as the largest index in tensor plus one. If not specified, no change of basis is applied.

Returns:

• ndarray –

  The dense tensor. Without a change of basis, its length along each axis is the largest index in tensor plus one; with a change of basis, it is the number of rows of basis. The entries have the type of the first value in tensor (before any change of basis).

Examples:

We convert a sparse tensor to a dense array. We then do the same in the basis consisting of the elements 1 and 3, which gives the same result as filter_tensor_indices with indices [1, 3].

>>> from cytools.utils import symmetric_sparse_to_dense
>>> tensor = {(0,1):0, (1,1):1, (1,2):2, (1,3):3, (2,3):4}
>>> symmetric_sparse_to_dense(tensor)
array([[0, 0, 0, 0],
       [0, 1, 2, 3],
       [0, 2, 0, 4],
       [0, 3, 4, 0]])
>>> basis = [[0,1,0,0],[0,0,0,1]]
>>> symmetric_sparse_to_dense(tensor, basis)
array([[1, 3],
       [3, 0]])

Source on GitHub

Convert a matrix given by its nonzero entries to a SciPy sparse matrix.

The input matrix can be given either as a list of triples [[a, b, M_ab], ...] or as a dictionary {(a, b): M_ab, ...}, where a and b are the row and column indices of an entry. The shape of the output is inferred from the largest row and column indices, and its entries are floats.

Parameters:

• arr (dict | list) –

  The matrix, either as a list of the form [[a, b, M_ab], ...] or as a dictionary of the form {(a, b): M_ab, ...}.

• sparse_type (str, default: 'dok' ) –

  The format of the output matrix. The options are "dok" (a scipy.sparse.dok_matrix) and "csr" (a scipy.sparse.csr_matrix).

Returns:

• dok_matrix | csr_matrix –

  The sparse matrix.

Raises:

• ValueError –

  If sparse_type is not one of the allowed options.

Examples:

We convert the 5x5 identity matrix into a sparse matrix, in the two available formats.

>>> from cytools.utils import to_sparse
>>> id_array = [[0,0,1],[1,1,1],[2,2,1],[3,3,1],[4,4,1]]
>>> M = to_sparse(id_array)
>>> M.format, M.shape, M.nnz
('dok', (5, 5), 5)
>>> to_sparse(id_array, sparse_type="csr").format
'csr'

The matrix can also be given as a dictionary.

>>> to_sparse({(0, 1): 2, (2, 0): 1}).toarray()
array([[0., 2.],
       [0., 0.],
       [1., 0.]])

Source on GitHub

Check whether the convex hull of the points is reflexive.

This is a minimal and fast check that avoids constructing a Polytope object. The polytope is considered reflexive if it is full-dimensional and all of its facets are at lattice distance 1 from the origin. Translations are not allowed, so the origin must be the unique interior point.

Parameters:

• points (ArrayLike) –

  The points whose convex hull defines the polytope.

• backend (str, default: 'qhull' ) –

  The backend used to compute the inequalities. The options are "ppl", "qhull", and "palp" (see poly_v_to_h).

Returns:

• bool –

  True if the convex hull of the points is reflexive, False otherwise.

Examples:

>>> from cytools.polytope import is_reflexive_barebones
>>> is_reflexive_barebones([[1,0],[0,1],[-1,-1]])
True
>>> is_reflexive_barebones([[2,1],[1,2],[0,0]])
False

Source on GitHub

Compute the H-representation of a polytope from its V-representation.

That is, map the points (or vertices) to the hyperplane inequalities defining their convex hull. Each row \((c_0, \ldots, c_{d-1}, c_d)\) of the output encodes the inequality

The points are expected to span a full-dimensional polytope. Otherwise, the output may contain equalities (PPL) or the computation may fail.

Parameters:

• pts (ArrayLike) –

  The input points, one per row.

• backend (str) –

  The backend to use. The options are "ppl", "qhull", and "palp".

Returns:

• tuple[ndarray, object] –

  The integer array of inequalities, one per row, and the convex hull object computed by the backend. The latter is a ppl.C_Polyhedron for PPL, a scipy.spatial.ConvexHull for QHull (or None for 1-dimensional polytopes), and None for PALP.

Raises:

• ValueError –

  If backend is not one of the allowed options.

Examples:

>>> from cytools.polytope import poly_v_to_h
>>> ineqs, _ = poly_v_to_h([[1,0],[0,1],[-1,-1]], backend="palp")
>>> ineqs
array([[ 2, -1,  1],
       [-1,  2,  1],
       [-1, -1,  1]])

Source on GitHub

Compute the lattice points of a polytope and the inequalities they saturate.

The lattice points contained in the convex hull of pts_in are enumerated within a bounding box, and for each of them the indices of the hyperplane inequalities that it saturates are found. The points are expected to span a full-dimensional polytope.

Parameters:

• pts_in (list[tuple]) –

  The points whose convex hull defines the polytope.

• ineqs (ArrayLike | None, default: None ) –

  The hyperplane inequalities defining the polytope, in the same format as the output of poly_v_to_h. If not specified, they are computed with backend.

• dim (int | None, default: None ) –

  The dimension of the polytope. If not specified, it is computed from the points.

• backend (str | None, default: None ) –

  The backend used to compute the inequalities if they are not given. The options are "ppl", "qhull", and "palp". If not specified, PPL is used for polytopes of dimension 1 to 4, and PALP otherwise. 0-dimensional polytopes always use PALP.

Returns:

• tuple[ndarray, list[frozenset[int]]] –

  The lattice points, one per row, and, for each point, the set of indices of the inequalities that it saturates. The points are not in any particular order.

Raises:

• ValueError –

  If the inequalities have coefficients that do not fit in a 32-bit integer.

Examples:

>>> from cytools.polytope import saturating_lattice_pts
>>> pts, sat = saturating_lattice_pts([(1,0),(0,1),(-1,-1)])
>>> pts
array([[-1, -1],
       [ 0,  0],
       [ 1,  0],
       [ 0,  1]], dtype=int32)
>>> sat
[frozenset({0, 1}), frozenset(), frozenset({1, 2}), frozenset({0, 2})]

Source on GitHub

Compute all triangulations of a point configuration using TOPCOM.

There is the option to only compute fine, regular, or star triangulations. All the triangulations are computed with TOPCOM (through triangulumancer) when the generator is first advanced, and they are then yielded one at a time. Regularity is checked for each triangulation as it is yielded.

Parameters:

• poly (Polytope) –

  The ambient polytope.

• pts (ArrayLike) –

  The points to be triangulated, specified by their labels in poly.

• only_fine (bool, default: False ) –

  Whether to restrict to fine triangulations.

• only_regular (bool, default: False ) –

  Whether to restrict to regular triangulations. Ignored if raw_output is true.

• only_star (bool, default: False ) –

  Whether to restrict to star triangulations.

• star_origin (int | None, default: None ) –

  The label of the point used as the star origin. It must be specified if only_star is true.

• backend (str | None, default: None ) –

  The optimizer used to check regularity. The available options are the backends of Cone.is_solid. If not specified, it is picked automatically.

• raw_output (bool, default: False ) –

  Whether to yield the triangulations as arrays of simplices (each specified by point labels) instead of as Triangulation objects. In this case regularity is not checked.

Yields:

• Triangulation | ndarray –

  The triangulations with the specified properties, or their arrays of simplices if raw_output is true.

Raises:

• ValueError –

  If pts is empty, or if only_star is true but star_origin is not specified.

Examples:

This function is not intended to be directly used, but it is used in the following example. We construct a polytope and find all of its triangulations. We try picking different restrictions and see how the number of triangulations changes.

>>> from cytools import Polytope
>>> p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],
...               [-1,-2,-1,-1],[-2,-1,-1,-1]])
>>> g = p.all_triangulations()
>>> next(g)  # takes some time while TOPCOM finds all the triangulations
A fine, regular, star triangulation of a 4-dimensional point...
>>> next(g)  # produces the next triangulation immediately
A fine, regular, star triangulation of a 4-dimensional point...
>>> len(p.all_triangulations(as_list=True))  # number of FRSTs
2
>>> len(p.all_triangulations(only_regular=False, only_star=False,
...                          only_fine=False, as_list=True))
6

Source on GitHub

Generate a fair sample of random regular triangulations.

The triangulations are fine and regular, and optionally star. This implements Algorithm #3 from the paper Bounding the Kreuzer-Skarke Landscape by Mehmet Demirtas, Liam McAllister, and Andres Rios-Tascon (arXiv:2008.01730).

This is a Markov chain Monte Carlo algorithm that involves taking random walks inside the subset of the secondary fan corresponding to fine triangulations and performing random flips. For details, please see Section 4.1 in the paper. The generator stops once N triangulations have been produced, once max_retries consecutive attempts fail to produce a new triangulation, or if no wall of the fine subset of the secondary fan can be found (in which case a message is printed).

Parameters:

• poly (Polytope) –

  The ambient polytope.

• pts (ArrayLike) –

  The points to be triangulated, specified by their labels in poly. They must be full-dimensional.

• N (int | None, default: None ) –

  Number of desired unique triangulations. If not specified, it generates as many triangulations as it can find until it has to retry max_retries times in a row to obtain a new triangulation.

• n_walk (int, default: 10 ) –

  Number of hit-and-run steps per triangulation.

• n_flip (int, default: 10 ) –

  Number of random flips performed per triangulation.

• initial_walk_steps (int, default: 20 ) –

  Number of hit-and-run steps to take before starting to record triangulations. Small values may result in a bias towards Delaunay-like triangulations.

• walk_step_size (float, default: 0.01 ) –

  Determines the size of random steps taken in the secondary fan. The algorithm may stall if it is too small.

• max_steps_to_wall (int, default: 10 ) –

  Maximum number of steps to take towards a wall of the subset of the secondary fan that corresponds to fine triangulations. If a wall is not found, a new random direction is selected. Setting this to be very large (>100) reduces performance. If this, or walk_step_size, is set too low, the algorithm may stall.

• fine_tune_steps (int, default: 8 ) –

  Number of steps to determine the location of a wall. Decreasing it improves performance, but might result in biased samples.

• max_retries (int, default: 50 ) –

  Maximum number of consecutive attempts to obtain a new triangulation (or to find a wall) before the process is terminated.

• make_star (bool, default: False ) –

  Whether to convert the obtained triangulations into star triangulations.

• backend (str, default: 'cgal' ) –

  The backend used to compute the triangulations. The available options are "cgal" and "qhull".

• seed (int | None, default: None ) –

  A seed for the random number generator. This can be used to obtain reproducible results.

Yields:

• Triangulation –

  The random triangulations with the specified properties.

Raises:

• Exception –

  If the point configuration is not full-dimensional.

Examples:

This function is not intended to be directly used, but it is used in the following example. We construct a polytope and find some random triangulations. The computation takes considerable time, but they should be a fair sample from the full set of triangulations (if the parameters are chosen correctly). For (some) machine learning purposes or when the fairness of the sample is not crucial, the fast sampler (Polytope.random_triangulations_fast) should be used instead.

>>> from cytools import Polytope
>>> p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],
...               [-1,-1,-1,-1]]).dual()
>>> g = p.random_triangulations_fair()
>>> next(g)  # takes a long time (a few minutes)
A fine, regular, star triangulation of a 4-dimensional point...
>>> next(g)  # also takes a few minutes
A fine, regular, star triangulation of a 4-dimensional point...
>>> rand_triangs = p.random_triangulations_fair(N=10, as_list=True)

It is worth noting that the time it takes to obtain each triangulation varies very significantly with the parameters used. The function tries to guess reasonable parameters, but it is better to adjust them to your desired balance between speed and fairness of the sampling.

Source on GitHub

Generate pseudorandom regular triangulations of a point configuration.

The triangulations are regular, and optionally fine and star. They are obtained by picking random heights around the Delaunay heights \(|p|^2\) from a Gaussian distribution. Duplicate triangulations are skipped. The generator stops once N triangulations have been produced, or once max_retries consecutive attempts fail to produce a new triangulation.

Parameters:

• poly (Polytope) –

  The ambient polytope.

• pts (ArrayLike) –

  The points to be triangulated, specified by their labels in poly.

• N (int | None, default: None ) –

  Number of desired unique triangulations. If not specified, it generates as many triangulations as it can find until it has to retry max_retries times in a row to obtain a new triangulation.

• c (float, default: 0.2 ) –

  The standard deviation of the Gaussian distribution used to pick the heights. A larger c results in a wider range of possible triangulations, but with a larger fraction of them being non-fine, which slows down the process when only_fine is true.

• max_retries (int, default: 500 ) –

  Maximum number of consecutive attempts to obtain a new triangulation before the process is terminated.

• make_star (bool, default: False ) –

  Whether to convert the obtained triangulations into star triangulations.

• only_fine (bool, default: True ) –

  Whether to restrict to fine triangulations.

• backend (str, default: 'cgal' ) –

  The backend used to compute the triangulations. The available options are "cgal" and "qhull".

• seed (int | None, default: None ) –

  A seed for the random number generator. This can be used to obtain reproducible results.

• verbosity (int, default: 0 ) –

  The verbosity level.

Yields:

• Triangulation –

  The random triangulations with the specified properties.

Examples:

This function is not intended to be directly used, but it is used in the following example. We construct a polytope and find some random triangulations. The triangulations are obtained very quickly, but they are not a fair sample of the space of triangulations.

>>> from cytools import Polytope
>>> p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],
...               [-1,-1,-1,-1]]).dual()
>>> g = p.random_triangulations_fast()
>>> next(g)  # runs very quickly
A fine, regular, star triangulation of a 4-dimensional point...
>>> next(g)  # keeps producing triangulations until it has trouble
A fine, regular, star triangulation of a 4-dimensional point...
>>> rand_triangs = p.random_triangulations_fast(N=10, as_list=True)

Source on GitHub

Convert between the hyperplane and ray descriptions of a cone.

The matrix M can be interpreted as the inward-pointing hyperplane normals of the cone \(\{x : Mx\geq 0\}\), in which case the returned matrix \(R\) contains rays that generate it, i.e. the cone is \(\{R^T\lambda : \lambda\geq 0\}\). Equivalently, by cone duality, M can be interpreted as the generating rays of the cone \(\{M^T\lambda : \lambda\geq 0\}\), in which case the returned matrix contains the hyperplane normals that define it. The computation is done with PPL, treating M as hyperplanes.

The output is a minimized description, but it is not guaranteed to consist of extremal rays. In particular, a lineality space is represented by pairs of opposite vectors \(\pm v\).

Parameters:

• M (ArrayLike) –

  The matrix defining the cone, one vector per row.

• verbosity (int, default: 0 ) –

  The verbosity level.

Returns:

• ndarray –

  The dual description of the cone, one vector per row. If there are no vectors, an empty array of shape (0, M.shape[1]) is returned.

Examples:

We find the rays of the cone defined by \(x_1\geq 0\) and \(x_1+x_2\geq 0\).

>>> from cytools.cone import dualize
>>> dualize([[1,0],[1,1]])
array([[ 0,  1],
       [ 1, -1]])

Source on GitHub

Solve the feasibility problem \(Ax\geq c\).

A point \(x\) is found such that \(a\cdot x\geq c\) for each row \(a\) of \(A\), minimizing the average of the \(a\cdot x\) (so that the point tends to lie close to the tip of the cone), and optionally subject to a lower bound on each coordinate.

Parameters:

• hyperplanes (ArrayLike) –

  The constraining hyperplanes, \(A\), one per row. The rows can also be sparse, given as objects whose items() method yields (index, value) pairs.

• c (float) –

  The "stretching", i.e. the minimum value of \(a\cdot x\) for each row \(a\).

• ambient_dim (int) –

  The ambient dimension, i.e. the number of columns of \(A\).

• backend (str) –

  The backend to use. The options are "highs" (an LP via HiGHS), "glop" (an LP via OR-Tools), "scip" (an integer program via OR-Tools), and "cpsat" (an integer program via OR-Tools CP-SAT).

• lower_bound (float | None, default: None ) –

  A lower bound on each of the variables. If not specified, the variables are unbounded below.

• verbose (bool, default: False ) –

  Whether to print extra diagnostic information.

Returns:

• ndarray | None –

  A feasible point, if one was found, or None otherwise. If there are no hyperplanes, a vector of ones is returned.

Examples:

We find a point satisfying \(x_1\geq 1\) and \(-x_1+x_2\geq 1\).

>>> from cytools.cone import feasibility
>>> feasibility([[1,0],[-1,1]], c=1, ambient_dim=2, backend="highs")
array([1., 2.])

Source on GitHub

Check whether a ray is extremal among a set of rays.

This is an auxiliary function used by Cone.extremal_rays. The ray R[i] is extremal if it cannot be written as a non-negative linear combination of the other rays. The output includes the index of the ray and any exception that was raised, so that it can be used when parallelizing.

Parameters:

• R (ArrayLike) –

  A matrix whose rows are the rays of the cone.

• i (int) –

  The index of the ray to check for extremality.

• extFlags (list[bool] | None, default: None ) –

  Flags indicating which rays are possibly extremal. Rays whose flag is false are known to be non-extremal and are excluded from the comparison. If not specified, all the other rays are used.

• method (str, default: 'lp' ) –

  The method used to check extremality. The options are "lp" (a linear program, recommended) and "nnls" (non-negative least squares).

• tol (float, default: 0.0001 ) –

  The tolerance for determining whether a ray is extremal. Only used if method is "nnls".

Returns:

• tuple[int, bool | None, Exception | None] –

  The index i, the truth value of the ray being extremal, and None. If an exception was raised during the computation, the truth value is None and the exception is returned as the third element.

Examples:

This function is not meant to be directly used by the end user. We check whether two rays are extremal.

>>> import numpy as np
>>> from cytools.cone import is_extremal
>>> R = np.array([[0,1],[1,1],[1,0]])
>>> is_extremal(R, 0)
(0, True, None)
>>> is_extremal(R, 1)
(1, False, None)

Source on GitHub

Enumerate the lattice points of a polytope.

The polytope is given by both its vertices and its inequalities, and need not be a lattice polytope. All lattice points in the bounding box of the vertices are tested against the inequalities, so this is only suitable for small polytopes. This is simpler than saturating_lattice_pts in the polytope module since it only enumerates the lattice points.

Parameters:

• verts (ArrayLike) –

  The vertices of the polytope, one per row. They may be rational.

• ineqs (ndarray) –

  The inequalities, as a matrix \(c\) where each row \(c_i\) encodes the inequality \(c_{i,0}x_0 + \dots + c_{i,d-1}x_{d-1} + c_{i,d} \geq 0\).

Returns:

• list[list[int]] –

  The lattice points of the polytope.

Examples:

>>> import numpy as np
>>> from cytools.h_polytope.h_polytope import lattice_points
>>> ineqs = np.array([[1,0,0],[0,1,0],[-6,-2,3]])
>>> lattice_points([[0,0],[0.5,0],[0,1.5]], ineqs)
[[0, 0], [0, 1]]

Source on GitHub

Compute the V-representation of a polytope from its H-representation.

That is, map hyperplane inequalities to the vertices of the polytope they define. The computation is done with PPL. If the inequalities are not integral, each row is first divided by the GCD of its entries and rounded, so rational inputs are made integral without changing the polytope.

Parameters:

• hypers (ArrayLike) –

  The inequalities, as a matrix \(c\) where each row \(c_i\) encodes the inequality \(c_{i,0}x_0 + \dots + c_{i,d-1}x_{d-1} + c_{i,d} \geq 0\).

• verbosity (int, default: 0 ) –

  The verbosity level.

Returns:

• tuple[ndarray, C_Polyhedron] –

  The vertices of the polytope, one per row, and the PPL polyhedron. The vertices have an integer dtype if they are all integral and a float dtype otherwise. If the inequalities have no solution, the array of vertices is empty.

Raises:

• ValueError –

  If the inequalities define an unbounded polyhedron.

Examples:

>>> from cytools.h_polytope.h_polytope import poly_h_to_v
>>> pts, poly = poly_h_to_v([[1,0,0],[0,1,0],[-6,-2,3]])
>>> pts
array([[0.5, 0. ],
       [0. , 1.5],
       [0. , 0. ]])

Source on GitHub

Check whether a newer version of CYTools is available.

The version number on the main branch of the CYTools GitHub repository is compared with the installed one. A message is printed if a newer version is available, and a warning is printed if the installed version is known to contain a serious bug. Nothing is printed otherwise.

This performs a network request, so it blocks until the request completes or times out (after 2 seconds). On import, CYTools instead runs it in a background daemon thread (see start_update_check), so that import cytools never waits on the network. Any failure (no network, DNS failure, timeout, ...) is silently ignored.

Examples:

We check for updates of CYTools. This is done automatically on import, so there is usually no need to do this.

>>> import cytools
>>> cytools.check_for_updates()

Source on GitHub

Start an update check in a background thread.

This runs check_for_updates in a background daemon thread and returns immediately. It is what runs on import, so that import cytools never blocks on the network. Since the thread is a daemon, it never keeps the interpreter alive, and any error it encounters is silently ignored.

The check is skipped entirely if the CYTOOLS_NO_UPDATE_CHECK environment variable is set to a non-empty value, which is useful on clusters and air-gapped machines.

Returns:

• Thread | None –

  The thread performing the check, or None if the check is disabled.

Examples:

We start a background update check. This is done automatically on import, so there is usually no need to do this.

>>> import cytools
>>> thread = cytools.start_update_check()

Source on GitHub