Skip to content

Cone

Cone

Cone(rays=None, hyperplanes=None, parse_inputs=True, check=True, copy=True, ambient_dim=None)

Rational polyhedral cone.

This class handles all computations relating to rational polyhedral cones, such as cone duality and extremal ray computations. It is mainly used for the study of Kähler and Mori cones. A cone can be specified either by a set of generating rays \(r_i\), in which case it is the set of non-negative linear combinations \(\sum_i \lambda_i r_i\) with \(\lambda_i\geq 0\), or by a set of inward-pointing hyperplane normals \(h_j\), in which case it is the set of points \(x\) with \(h_j\cdot x\geq 0\) for all \(j\). The other description is computed on demand.

The input rays or hyperplane normals are stored as integer vectors. Floating-point input is converted by dividing each row by the (numerical) GCD of its entries and rounding. This class is also available at the root of the package, so it can be imported with from cytools import Cone.

Parameters:

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

    The rays that generate the cone, one per row. If it is not specified, then hyperplanes must be specified.

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

    The inward-pointing normals to the hyperplanes that define the cone, one per row. If it is not specified, then rays must be specified. An empty list defines the whole ambient space.

  • parse_inputs (bool, default: True ) –

    Whether to parse and normalize the input. If false, hyperplanes must be a 2D NumPy array of integers, and it is stored as is. This is only supported when constructing a cone from hyperplanes.

  • check (bool, default: True ) –

    Whether to divide each integer input row by the GCD of its entries, dropping rows of zeros (with a warning). Floating-point input is always normalized this way. Recommended when constructing a cone directly.

  • copy (bool, default: True ) –

    Whether to copy the input data. If false, the input array may be used directly when no conversion is needed. Recommended.

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

    The dimension of the ambient space. It is only used when hyperplanes is empty, in which case it must be specified unless it can be inferred from the shape of the input.

Raises:

  • ValueError –

    If not exactly one of rays or hyperplanes is specified, if the input is not a 2D matrix, if the ambient dimension is zero, or if hyperplanes is empty and ambient_dim is missing or inconsistent with its shape.

  • NotImplementedError –

    If parse_inputs is false and rays are given, or if the input has an unsupported data type.

  • Exception –

    If arbitrary-precision (fmpz or fmpq) input is given without enabling the experimental features.

Warning

This class is primarily tailored to pointed (i.e. strongly convex) cones. There are a few computations, such as finding extremal rays, that may produce some unexpected results when working with non-pointed cones.

Examples:

We construct a cone in two different ways, first from a list of rays and then from a list of hyperplane normals, and verify that the two inputs result in the same cone.

>>> from cytools import Cone
>>> c1 = Cone([[0,1],[1,1]])  # equivalently, Cone(rays=[[0,1],[1,1]])
>>> c2 = Cone(hyperplanes=[[1,0],[-1,1]])
>>> c1 == c2
True

Methods:

  • ambient_dimension –

    Return the dimension of the ambient lattice.

  • clear_cache –

    Clear the cached results of any previous computation.

  • contains –

    Check whether points, or another cone, are contained in the cone.

  • dimension –

    Return the dimension of the cone.

  • dual_cone –

    Return the dual cone.

  • extremal_hyperplanes –

    Return the extremal hyperplanes of the cone.

  • extremal_rays –

    Return the extremal rays of the cone.

  • faces –

    Return the faces of the cone, organized by codimension.

  • facets –

    Return the facets of the cone.

  • find_grading_vector –

    Find a grading vector for the cone.

  • find_interior_point –

    Find a point in the strict interior of the cone.

  • find_lattice_points –

    Find lattice points in the cone.

  • hilbert_basis –

    Return the Hilbert basis of the cone.

  • hyperplanes –

    Return the inward-pointing normals to the hyperplanes that define the cone.

  • intersection –

    Compute the intersection with another cone, or with a list of cones.

  • is_degenerate –

    Check whether the cone is degenerate.

  • is_pointed –

    Return whether the cone is pointed (i.e. strongly convex).

  • is_simplicial –

    Return whether the cone is simplicial.

  • is_smooth –

    Return whether the cone is smooth.

  • is_solid –

    Return whether the cone is solid, i.e. whether it is full-dimensional.

  • lineality_space –

    Return the lineality space of the cone as a Cone object.

  • pointed_space –

    Return the pointed part of the cone.

  • rays –

    Return the (not necessarily extremal) rays that generate the cone.

  • tip_of_stretched_cone –

    Find the tip of the stretched cone.

ambient_dimension

ambient_dimension()

Return the dimension of the ambient lattice.

Also available as ambient_dim.

Returns:

  • int –

    The dimension of the ambient lattice.

Examples:

We construct a cone and find the dimension of the ambient lattice.

>>> c = Cone([[0,1,0],[1,1,0]])
>>> c.ambient_dimension()
3

clear_cache

clear_cache()

Clear the cached results of any previous computation.

The input data (rays or hyperplane normals, whichever was used to construct the cone) is kept, while everything that was derived from it, including the other description of the cone, is discarded.

Examples:

We construct a cone, compute its extremal rays, clear the cache and then compute them again.

>>> c = Cone([[1,0],[1,1],[0,1]])
>>> c.extremal_rays()
array([[1, 0],
       [0, 1]])
>>> c.clear_cache()  # clears the cached result
>>> c.extremal_rays()  # the extremal rays are recomputed
array([[1, 0],
       [0, 1]])

contains

contains(other, eps=0)

Check whether points, or another cone, are contained in the cone.

A point \(x\) is considered to be contained if \(h\cdot x\geq \epsilon\) for every hyperplane normal \(h\) returned by hyperplanes. With the default \(\epsilon=0\) this checks membership in the (closed) cone, while a positive \(\epsilon\) can be used to check whether a point is in the strict interior. A cone is contained if all of its generating rays are.

Parameters:

  • other (ArrayLike | Cone) –

    The object to check. A 1D array is treated as a single point, a 2D array as a list of points (one per row), and a Cone is checked through its rays.

  • eps (float, default: 0 ) –

    The minimum value that \(h\cdot x\) must attain for every hyperplane normal \(h\).

Returns:

  • bool | tuple[bool, ...] –

    Whether the point or the cone is contained. If a 2D array is given, a tuple with one truth value per point is returned.

Examples:

We construct a cone and check whether some points and a cone are contained in it.

>>> c = Cone([[1,0],[0,1]])
>>> bool(c.contains([1,1]))
True
>>> [bool(b) for b in c.contains([[1,1],[-1,0]])]
[True, False]
>>> bool(c.contains([1,0], eps=1))  # on the boundary
False
>>> c.contains(Cone([[1,1],[1,2]]))
True

dimension

dimension()

Return the dimension of the cone.

This is the dimension of the linear span of the cone. If the cone was constructed from hyperplanes and is not solid, its rays are computed first, which can be expensive. The result is cached.

Also available as dim.

Returns:

  • int –

    The dimension of the cone.

Examples:

We construct a cone and find its dimension.

>>> c = Cone([[0,1,0],[1,1,0]])
>>> c.dimension()
2

dual_cone

dual_cone()

Return the dual cone.

The dual cone is the set of points \(y\) such that \(y\cdot x\geq 0\) for all points \(x\) in the cone. It is constructed by using the rays of this cone as the hyperplane normals of the dual (or vice versa), and it is cached. The dual of the returned cone is this cone.

Also available as dual.

Returns:

  • Cone –

    The dual cone.

Examples:

We construct a cone and find its dual cone.

>>> c = Cone([[0,1],[1,1]])
>>> c.dual_cone()
A rational polyhedral cone in RR^2 defined by 2 hyperplanes
>>> c.dual_cone().rays()
array([[ 1,  0],
       [-1,  1]])

extremal_hyperplanes

extremal_hyperplanes(tol=0.0001, minimal=True, method='extremalrays', verbose=False)

Return the extremal hyperplanes of the cone.

These are the extremal rays of the dual cone, i.e. a set of inward-pointing hyperplane normals that define the cone and are not redundant. This is computed with extremal_rays on the dual cone, so see its documentation for details.

Parameters:

  • tol (float, default: 0.0001 ) –

    The tolerance for deciding whether a hyperplane is extremal. Only used if method is "nnls".

  • minimal (bool, default: True ) –

    Whether to return a minimal set of hyperplanes. This only matters when the dual cone is not pointed, i.e. when the cone is not solid.

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

    The backend used to prune the hyperplanes. The options are "extremalrays", "legacy" (or its synonym "lp"), and "nnls".

  • verbose (bool, default: False ) –

    Whether to print the progress of the computation.

Returns:

  • ndarray –

    The extremal hyperplane normals of the cone, one per row.

Examples:

We construct a cone and find its extremal hyperplanes.

>>> c = Cone([[0,1],[1,1],[1,0]])
>>> c.extremal_hyperplanes()
array([[0, 1],
       [1, 0]])

extremal_rays

extremal_rays(tol=0.0001, minimal=True, method='extremalrays', verbose=False)

Return the extremal rays of the cone.

A ray is extremal if it cannot be written as a non-negative linear combination of the other rays. For pointed cones, the extremal rays form the unique minimal generating set of the cone. For non-pointed cones there is no such unique set, and the behaviour depends on minimal.

Duplicate rays are removed first, and the order of the returned rays need not match the order of the input rays. The result is cached separately for each value of minimal (the cache does not depend on method or tol). If the cone was constructed from hyperplanes, the extremal rays are also stored as its generating rays.

Parameters:

  • tol (float, default: 0.0001 ) –

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

  • minimal (bool, default: True ) –

    Whether to return a minimal generating set of rays. This only matters for non-pointed cones. If true, the cone is decomposed into its lineality space and its pointed part (see lineality_space and pointed_space), and the extremal rays of both are returned together, so the result need not be a subset of the generating rays. If false, the generating rays that are not non-negative combinations of the others are returned, which need not form a minimal generating set.

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

    The backend used to prune the rays. The options are: "extremalrays", which uses the extremalrays package and falls back to "legacy" if the package is not installed (with a warning) or if the cone is not pointed; "legacy" (or its synonym "lp"), which checks each ray by solving a linear program; and "nnls", which checks each ray by solving a non-negative least-squares problem.

  • verbose (bool, default: False ) –

    Whether to print the progress of the computation.

Returns:

  • ndarray –

    The extremal rays of the cone, one per row.

Raises:

  • ValueError –

    If method is not one of the allowed options.

  • RuntimeError –

    If the check of a ray fails repeatedly with the "legacy" or "nnls" backends.

Note

The "legacy" and "nnls" backends check the rays in parallel for pointed cones with at least 32 rays, using as many CPU threads as are available. To fix the number of threads, you can set the n_threads variable in the cytools.config module.

Examples:

We construct a cone and find its extremal rays.

>>> c = Cone([[0,1],[1,1],[1,0]])
>>> c.extremal_rays()
array([[0, 1],
       [1, 0]])

faces

faces(codim=None, include_self=False, verbosity=0)

Return the faces of the cone, organized by codimension.

The faces are found from the incidences between the extremal rays and the extremal hyperplanes of the cone. Each face of positive dimension is returned as a Cone generated by its extremal rays, while the zero-dimensional face (the origin) is returned as a Cone defined by hyperplanes. The result is cached.

Only pointed cones are currently supported. This computation can be expensive for high-dimensional cones, and it is slower than facets if only the facets are needed.

Parameters:

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

    The codimension of the desired faces. If not specified, the faces of all codimensions are returned. When set to 0, the cone itself is returned.

  • include_self (bool, default: False ) –

    Whether to include the cone itself (the codimension-0 face) when returning the faces of all codimensions.

  • verbosity (int, default: 0 ) –

    The verbosity level.

Returns:

  • tuple[Cone, ...] | tuple[tuple[Cone, ...], ...] –

    If codim is specified, the faces of that codimension. Otherwise, a tuple with one tuple of faces per codimension, in increasing order, from codimension 1 (or 0 if include_self is true) up to the dimension of the cone.

Raises:

Examples:

We construct a simplicial cone and find the number of faces of each codimension, as well as the rays of one of its facets.

>>> c = Cone([[1,0,0],[0,1,0],[0,0,1]])
>>> [len(f) for f in c.faces()]
[3, 3, 1]
>>> c.faces(codim=1)[0].rays()
array([[1, 0, 0],
       [0, 0, 1]])

facets

facets(verbosity=0)

Return the facets of the cone.

The facets are the faces of codimension 1. If the cone is simplicial, they are the cones generated by each subset of all but one of the extremal rays. Otherwise, they are computed by saturating, one at a time, each extremal hyperplane that does not vanish on the whole cone, which requires computing the extremal hyperplanes. In this case, the facets are returned as cones defined by hyperplanes. The result is not cached.

Parameters:

  • verbosity (int, default: 0 ) –

    The verbosity level.

Returns:

  • list[Cone] –

    The facets of the cone.

Examples:

We construct a non-simplicial cone and find its facets.

>>> c = Cone([[1,0,1],[0,1,1],[-1,0,1],[0,-1,1]])
>>> f = c.facets()
>>> len(f)
4
>>> f[0].extremal_rays()
array([[-1,  0,  1],
       [ 0,  1,  1]])

find_grading_vector

find_grading_vector(backend=None)

Find a grading vector for the cone.

A grading vector is a vector \(v\) such that \(v\cdot p>0\) for every non-zero point \(p\) in the cone. Thus, it must be in the strict interior of the dual cone, so it only exists for pointed cones. This function returns an integral grading vector, found with find_interior_point on the dual cone.

Parameters:

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

    The optimizer to use. The options are the same as for find_interior_point. If not specified, the default of that function is used.

Returns:

  • ndarray | None –

    An integral grading vector. If it could not be found, then None is returned.

Raises:

Examples:

We construct a cone and find a grading vector.

>>> c = Cone([[3,2],[5,3]])
>>> c.find_grading_vector()
array([-1,  2])

find_interior_point

find_interior_point(c=1, lower=None, integral=False, backend=None, check=True, show_hints=False, verbose=False)

Find a point in the strict interior of the cone.

If the rays of the cone are known and neither backend nor lower is specified, this simply returns the sum of the rays divided by their GCD (and, unless integral is true, divided by the number of rays), or None if the cone is not solid. In this case, c and check are ignored.

Otherwise, an optimization problem is solved to find a point \(x\) with \(h\cdot x\geq c\) for every hyperplane normal \(h\) returned by hyperplanes. The LP and integer programming backends minimize the average of the \(h\cdot x\) (see feasibility), while the quadratic programming backends find the tip of the stretched cone (see tip_of_stretched_cone).

Parameters:

  • c (float, default: 1 ) –

    A positive number specifying the stretching of the cone, i.e. the minimum value of \(h\cdot x\) for each hyperplane normal \(h\). Not used if the point is obtained from the rays.

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

    A lower bound on each of the coordinates of the point. It is not supported by the "mosek", "osqp", and "cvxopt" backends.

  • integral (bool, default: False ) –

    Whether the point should have integral coordinates. If the optimizer returns a non-integral point, its multiples are rounded until one that is strictly interior is found.

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

    The optimizer to use. The options are "highs" (an LP via HiGHS), "glop" (an LP via OR-Tools), "scip" (an integer program via OR-Tools), "cpsat" (an integer program via OR-Tools CP-SAT), and "mosek", "osqp", and "cvxopt" (quadratic programs). If not specified and the rays are not known, "mosek" is used if it is activated and the ambient dimension is at least 25, and "highs" otherwise.

  • check (bool, default: True ) –

    Whether to verify that the point found by the optimizer is in the strict interior of the cone.

  • show_hints (bool, default: False ) –

    Whether to print hints about possible reasons for failure of the quadratic programming backends.

  • verbose (bool, default: False ) –

    Whether to print extra diagnostic information.

Returns:

  • ndarray | None –

    A point in the strict interior of the cone. If no point is found, then None is returned.

Raises:

  • ValueError –

    If an invalid backend is specified, or if lower is specified with a backend that does not support it.

Examples:

We construct a cone and find some interior points.

>>> c = Cone([[3,2],[5,3]])
>>> c.find_interior_point()
array([4. , 2.5])
>>> c.find_interior_point(integral=True)
array([8, 5])

find_lattice_points

find_lattice_points(min_points=None, max_deg=None, grading_vector=None, c=0, max_coord=None, deg_window=0, filter_function=None, process_function=None, fast_mode=True, max_B=10000, verbose=False)

Find lattice points in the cone.

The lattice points \(x\) that are found satisfy \(h\cdot x\geq c\) for every hyperplane normal \(h\) returned by hyperplanes, so with the default c=0 they are the lattice points of the cone. Either the minimum number of points to find or the maximum degree of the points must be specified. There are two algorithms:

  • A fast mode, which is used when fast_mode is true, min_points is specified, the ambient dimension is at most 10, and none of max_deg, grading_vector, filter_function, or process_function are given. It enumerates the lattice points in boxes of increasing size (i.e. bounded by the \(L^\infty\) norm) using the latticepts package, and it returns at least min_points points, unless the box size reaches max_B. In this mode, max_coord and deg_window are ignored, and c is converted to integers.
  • The general algorithm, which requires the cone to be pointed. It uses OR-Tools CP-SAT to find the points whose degree, i.e. their dot product with the grading vector, is at most max_deg. If max_deg is not given, it instead searches in windows of deg_window + 1 consecutive degrees, starting from degree 0, until at least min_points points are found. All the points in the last window are kept, so more than min_points points may be returned. If both max_deg and min_points are given, min_points is ignored. The points are returned sorted by degree, and lexicographically within each degree.

Parameters:

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

    The minimum number of points to find. Note that this count does not take filter_function into account.

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

    The maximum degree of the points to find. This is useful when working with a preferred grading vector.

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

    The grading vector that will be used. If not specified, then it is computed with find_grading_vector.

  • c (float | ArrayLike, default: 0 ) –

    The minimum allowed stretching. It can be a single number, or one number per hyperplane, in the order of hyperplanes.

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

    The maximum absolute value of the coordinates of the points in the general algorithm. If not specified, the largest bound supported by CP-SAT is used.

  • deg_window (int, default: 0 ) –

    The number of additional degrees that are searched at a time when using min_points in the general algorithm. The searched windows are the degrees from \(n(w+1)\) to \(n(w+1)+w\), for \(n=0,1,2,\dots\) and \(w\) equal to deg_window.

  • filter_function (Callable | None, default: None ) –

    A function that takes a point (as a tuple of integers) and returns a boolean indicating whether to keep it.

  • process_function (Callable | None, default: None ) –

    A function that is called on each point (as a tuple of integers) as it is found. This is useful to avoid first constructing a large list of points and then processing it. When it is specified, nothing is returned.

  • fast_mode (bool, default: True ) –

    Whether to allow the fast mode described above.

  • max_B (int, default: 10000 ) –

    The maximum box size (i.e. \(L^\infty\) norm) allowed in the fast mode.

  • verbose (bool, default: False ) –

    Whether to print extra diagnostic information.

Returns:

  • ndarray | None –

    The lattice points, one per row. None is returned if process_function is specified, or if the search with max_deg was unsuccessful (in which case a message is printed).

Raises:

  • Exception –

    If neither min_points nor max_deg is specified, if both filter_function and process_function are specified, if the general algorithm is used on a cone that is not pointed, or if the grading vector is not valid.

  • ValueError –

    If max_coord is specified and not enough points can be found within the corresponding box.

Examples:

We construct a cone and find at least 20 lattice points in it.

>>> c = Cone([[3,2],[5,3]])
>>> pts = c.find_lattice_points(min_points=20)
>>> len(pts)
20

We can also find all the points up to a given degree. Here the grading vector is computed automatically.

>>> c.find_lattice_points(max_deg=2)
array([[ 0,  0],
       [ 3,  2],
       [ 5,  3],
       [ 6,  4],
       [ 8,  5],
       [10,  6]])

Let's also give an example where we use a function to apply some filtering. This can be something very complicated, but here we just pick the points where all coordinates are odd.

>>> def odd_coords(pt):
...     return all(c%2 for c in pt)
>>> pts = c.find_lattice_points(min_points=20, filter_function=odd_coords)
>>> len(pts)  # now we get only 6 points
6

Finally, let's give an example where we process the data as it comes instead of first constructing a list. In this simple example we just print each point with odd coordinates, but in general it can be a complex algorithm.

>>> def process_function(pt):
...     if all(c%2 for c in pt):
...         print(f"Processing point {pt}")
>>> c.find_lattice_points(min_points=20, process_function=process_function)
Processing point (5, 3)
Processing point (11, 7)
Processing point (15, 9)
Processing point (17, 11)
Processing point (21, 13)
Processing point (25, 15)

hilbert_basis

hilbert_basis()

Return the Hilbert basis of the cone.

The Hilbert basis is the minimal set of lattice points that generate all the lattice points in the cone as non-negative integer combinations. Normaliz is used for the computation, so its normaliz executable must be installed and available on the PATH. The result is cached.

Returns:

  • ndarray –

    The vectors forming the Hilbert basis, one per row.

Raises:

Examples:

We compute the Hilbert basis of a two-dimensional cone.

>>> c = Cone([[1,3],[2,1]])
>>> c.hilbert_basis()
array([[1, 1],
       [1, 2],
       [1, 3],
       [2, 1]])

hyperplanes

hyperplanes(use_extremal_rays=False, verbosity=0)

Return the inward-pointing normals to the hyperplanes that define the cone.

The cone is the set of points \(x\) such that \(h\cdot x\geq 0\) for every returned normal \(h\). If the cone was constructed from hyperplane normals, these are returned (after the normalization done at construction). Otherwise, they are computed from the rays with PPL and cached. In that case, if the cone is not solid, the linear subspace it spans is cut out by pairs of opposite normals \(\pm h\). This can be slow for cones in more than about 12 dimensions.

Parameters:

  • use_extremal_rays (bool, default: False ) –

    Whether to first compute the extremal rays and use them in this computation, instead of all the generating rays. Only relevant if the hyperplanes are not already known.

  • verbosity (int, default: 0 ) –

    The verbosity level.

Returns:

  • ndarray –

    A copy of the inward-pointing hyperplane normals, one per row.

Examples:

We construct two cones and find their hyperplane normals.

>>> c1 = Cone([[0,1],[1,1]])
>>> c2 = Cone(hyperplanes=[[0,1],[1,1]])
>>> c1.hyperplanes()
array([[ 1,  0],
       [-1,  1]])
>>> c2.hyperplanes()
array([[0, 1],
       [1, 1]])

intersection

intersection(other)

Compute the intersection with another cone, or with a list of cones.

The resulting cone is defined by the combined hyperplane normals of all the cones.

Parameters:

  • other (Cone | Iterable[Cone]) –

    The other cone that is being intersected, or a list of cones to intersect with.

Returns:

  • Cone –

    The cone that results from the intersection.

Raises:

  • ValueError –

    If other is a list containing objects that are not cones, or cones with a different ambient dimension.

Examples:

We construct two cones and find their intersection.

>>> c1 = Cone([[1,0],[1,2]])
>>> c2 = Cone([[0,1],[2,1]])
>>> c3 = c1.intersection(c2)
>>> c3.rays()
array([[2, 1],
       [1, 2]])

is_degenerate

is_degenerate(use_extremal_hyperplanes=True, M=None, certificate=False, verbosity=0)

Check whether the cone is degenerate.

For a cone \(\{x : Hx\geq 0\}\) in \(d\) ambient dimensions, this checks whether some non-zero \(x\) in the cone saturates (i.e. satisfies \(h\cdot x=0\)) at least \(d+1\) of the hyperplanes \(h\) simultaneously. If so, the cone is considered degenerate. It is more difficult to compute the rays of a degenerate cone.

This is representation-sensitive: a cone may be degenerate for a certain matrix of hyperplane normals \(H\) but not for another one that defines the same cone. It is therefore best to use the extremal hyperplanes.

The check first tries \(x=(1,\ldots,1)\), and otherwise searches for such an \(x\) among the integer vectors with OR-Tools CP-SAT.

Parameters:

  • use_extremal_hyperplanes (bool, default: True ) –

    Whether to check the extremal hyperplanes for degeneracy. If false, the hyperplanes returned by hyperplanes are used.

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

    The bound on the absolute value of the coordinates of \(x\) in the search. If not specified, the largest bound supported by CP-SAT is used.

  • certificate (bool, default: False ) –

    Whether to also return the point \(x\) that was found and the hyperplanes that the solver claims it saturates.

  • verbosity (int, default: 0 ) –

    The verbosity level.

Returns:

  • bool | tuple | None –

    Whether the cone is degenerate, or None if the solver failed. If certificate is true, a tuple (degenerate, (x, z)) is returned instead, where x is the point that was found and z is an array of 0s and 1s indicating which hyperplanes it saturates (both can be None if no point was found).

Examples:

We construct a cone over a square and a cone over a pyramid with a pentagonal base, and check whether they are degenerate. The apex of the pyramid lies on five facets.

>>> c1 = Cone([[1,0,1],[0,1,1],[-1,0,1],[0,-1,1]])
>>> bool(c1.is_degenerate())
False
>>> c2 = Cone([[1,0,0,1],[0,1,0,1],[-1,1,0,1],[-1,0,0,1],[0,-1,0,1],
...            [0,0,1,1]])
>>> bool(c2.is_degenerate())
True

is_pointed

is_pointed(backend='dual', tol=1e-07)

Return whether the cone is pointed (i.e. strongly convex).

A cone is pointed if it contains no line, i.e. if there is no non-zero \(x\) such that both \(x\) and \(-x\) are in the cone. The available backends are:

  • "dual": checks whether the dual cone is solid (see is_solid). This works for any cone.
  • "null": checks whether the matrix of hyperplane normals \(H\) has full rank, since otherwise any non-zero \(x\) with \(Hx=0\) spans a line in the cone. This requires the hyperplanes to be known.
  • "lp" and "nnls": check whether there are non-negative coefficients \(\lambda\neq 0\) such that \(\sum_i\lambda_i r_i=0\) for the generating rays \(r_i\), using a linear program or non-negative least squares, respectively. These can only be used when the hyperplanes are not known.

The result is cached, regardless of the backend used.

Also available as is_strongly_convex.

Parameters:

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

    The backend to use. The options are "dual", "null", "lp", and "nnls", as described above.

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

    The tolerance for determining when a linear subspace is found. This is only used by the "nnls" backend.

Returns:

  • bool –

    The truth value of the cone being pointed.

Raises:

  • ValueError –

    If the backend is unknown, or if it cannot be used with the known description of the cone.

Examples:

We construct two cones and check if they are pointed.

>>> c1 = Cone([[1,0],[0,1]])
>>> c2 = Cone([[1,0],[0,1],[-1,0]])
>>> c1.is_pointed()
True
>>> c2.is_pointed()
False

is_simplicial

is_simplicial()

Return whether the cone is simplicial.

A cone is simplicial if its extremal rays are linearly independent, i.e. if the number of extremal rays equals the dimension of the cone. If the rays are not known and the cone is solid, it is instead checked whether the number of extremal hyperplanes equals the dimension, which is equivalent in that case. The result is cached.

Returns:

  • bool –

    The truth value of the cone being simplicial.

Examples:

We construct two cones and check if they are simplicial.

>>> c1 = Cone([[1,0,0],[0,1,0],[0,0,1]])
>>> c2 = Cone([[1,0,0],[0,1,0],[0,0,1],[1,1,-1]])
>>> c1.is_simplicial()
True
>>> c2.is_simplicial()
False

is_smooth

is_smooth()

Return whether the cone is smooth.

A cone is smooth if it is simplicial and its extremal rays either form a basis of the ambient lattice, or they can be extended into one. The result is cached.

Returns:

  • bool –

    The truth value of the cone being smooth.

Examples:

We construct two cones and check if they are smooth.

>>> c1 = Cone([[1,0,0],[0,1,0],[0,0,1]])
>>> c2 = Cone([[2,0,1],[0,1,0],[1,0,2]])
>>> c1.is_smooth()
np.True_
>>> c2.is_smooth()
np.False_

is_solid

is_solid(backend=None)

Return whether the cone is solid, i.e. whether it is full-dimensional.

If the generating rays are known, this is simply checked by computing the dimension of the linear space that they span, and backend is ignored. When only the hyperplane inequalities are known this can be a more difficult problem. When using PPL as the backend, the cone is explicitly constructed and its dimension is checked. The other backends try to find a point in the strict interior of the cone with find_interior_point, which fails if the cone is not solid. The latter approach is much faster, but there could be extremely narrow cones where the optimization fails and this function returns a false negative. The result is cached.

Also available as is_full_dimensional.

Parameters:

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

    The backend to use when the rays are not known. The options are "ppl", and any of the backends of find_interior_point. If not specified, the default backend of find_interior_point is used.

Returns:

  • bool –

    The truth value of the cone being solid.

Raises:

  • ValueError –

    If an invalid backend is specified.

Examples:

We construct two cones and check if they are solid.

>>> c1 = Cone([[1,0],[0,1]])
>>> c2 = Cone([[1,0,0],[0,1,0]])
>>> c1.is_solid()
True
>>> c2.is_solid()
False

lineality_space

lineality_space()

Return the lineality space of the cone as a Cone object.

The lineality space is the largest linear subspace contained in the cone, i.e. the set of \(x\) such that \(Hx=0\), where \(H\) is the matrix of hyperplane normals. Since it is a linear subspace, packaging it as a cone is somewhat formal: it is returned as the cone generated by an integral basis of the subspace together with minus the sum of the basis vectors, and these rays are also set as its extremal rays. For a pointed cone, the lineality space is the origin.

Returns:

  • Cone –

    The lineality space of the cone.

Examples:

We construct a cone that is a half-plane and find its lineality space.

>>> c = Cone([[1,0],[0,1],[-1,0]])
>>> c.lineality_space().rays()
array([[ 1,  0],
       [-1,  0]])

pointed_space

pointed_space()

Return the pointed part of the cone.

A cone can be decomposed as the sum of its lineality space (see lineality_space) and a pointed cone. The pointed part is obtained by intersecting the cone with the orthogonal complement of the lineality space, and it is returned as a cone defined by hyperplanes.

Returns:

  • Cone –

    The pointed part of the cone.

Examples:

We construct a cone that is a half-plane and find its pointed part.

>>> c = Cone([[1,0],[0,1],[-1,0]])
>>> c.pointed_space().rays()
array([[0, 1]])

rays

rays(use_extremal_hyperplanes=False, verbosity=0)

Return the (not necessarily extremal) rays that generate the cone.

If the cone was constructed from rays, these are returned (after the normalization done at construction). Otherwise, they are computed from the hyperplane normals with PPL and cached. In that case, a lineality space is represented by pairs of opposite rays \(\pm v\). This can be slow for cones in more than about 12 dimensions.

Parameters:

  • use_extremal_hyperplanes (bool, default: False ) –

    Whether to first compute the extremal hyperplanes and use them in this computation, instead of all the defining hyperplanes. Only relevant if the rays are not already known.

  • verbosity (int, default: 0 ) –

    The verbosity level.

Returns:

  • ndarray –

    A copy of the rays that generate the cone, one per row.

Examples:

We construct two cones and find their generating rays.

>>> c1 = Cone([[0,1],[1,1]])
>>> c2 = Cone(hyperplanes=[[0,1],[1,1]])
>>> c1.rays()
array([[0, 1],
       [1, 1]])
>>> c2.rays()
array([[ 1,  0],
       [-1,  1]])

tip_of_stretched_cone

tip_of_stretched_cone(c=1, backend=None, check=True, constraint_error_tol=0.05, max_iter=10 ** 6, show_hints=True, verbose=False)

Find the tip of the stretched cone.

The stretched cone is the region of points \(x\) in the cone such that \(h\cdot x\geq c\) for every hyperplane normal \(h\) returned by hyperplanes. Its tip is the point in this region with the smallest (Euclidean) norm. Note that the normals are integral vectors that are not normalized, so \(c\) is not literally a distance to the hyperplanes.

If the cone has no defining hyperplanes (i.e. it is the whole space), a vector of ones is returned.

Parameters:

  • c (float, default: 1 ) –

    A positive number specifying the stretching of the cone.

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

    The optimizer to use. The options are "mosek", "osqp", and "cvxopt", which solve the quadratic program, and "highs" and "glop", which solve a linear program whose solution is only an approximation of the tip (see find_interior_point). If not specified, "osqp" is used for ambient dimension \(d<25\). For \(d\geq 25\), "mosek" is used if it is activated, and "highs" otherwise.

  • check (bool, default: True ) –

    Whether to check that the output of the optimizer satisfies the constraints within constraint_error_tol.

  • constraint_error_tol (float, default: 0.05 ) –

    The error tolerance for the linear constraints.

  • max_iter (int, default: 10 ** 6 ) –

    The maximum number of iterations allowed for the quadratic programming backends. If this function is returning None, then increasing this parameter (up to \(2^{31}-1\)) might resolve the issue. It is ignored by the "highs" and "glop" backends.

  • show_hints (bool, default: True ) –

    Whether to print hints about possible reasons for failure when no solution is found.

  • verbose (bool, default: False ) –

    Whether to print extra diagnostic information.

Returns:

  • ndarray | None –

    The location of the tip. If it could not be found, or if it does not pass the check, then None is returned.

Raises:

  • ValueError –

    If an invalid backend is specified.

  • Exception –

    If "mosek" is specified but Mosek is not activated.

Note

Finding the tip requires quadratic programming, since the norm of a vector is being minimized. For dimensions up to around 25, this can easily be done with open-source solvers like OSQP or CVXOPT. For higher dimensions this becomes a difficult task that only the Mosek optimizer handles well. However, Mosek is closed-source and requires a license. For this reason, when Mosek is not activated, the default for high dimensions is to solve a linear problem with HiGHS, which finds a good approximation of the tip.

Examples:

We construct two cones and find the locations of the tips of the stretched cones.

>>> c1 = Cone([[1,0],[0,1]])
>>> c2 = Cone([[3,2],[5,3]])
>>> c1.tip_of_stretched_cone(1)
array([1., 1.])
>>> c2.tip_of_stretched_cone(1)
array([8., 5.])