Module bayesvalidrox.surrogate_models.glexindex
Multi indices for monomial exponents. Credit: Jonathan Feinberg https://github.com/jonathf/numpoly/blob/master/numpoly/utils/glexindex.py
Expand source code
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Multi indices for monomial exponents.
Credit: Jonathan Feinberg
https://github.com/jonathf/numpoly/blob/master/numpoly/utils/glexindex.py
"""
import numpy
import numpy.typing
def glexindex(start, stop=None, dimensions=1, cross_truncation=1.,
graded=False, reverse=False):
"""
Generate graded lexicographical multi-indices for the monomial exponents.
Args:
start (Union[int, numpy.ndarray]):
The lower order of the indices. If array of int, counts as lower
bound for each axis.
stop (Union[int, numpy.ndarray, None]):
The maximum shape included. If omitted: stop <- start; start <- 0
If int is provided, set as largest total order. If array of int,
set as upper bound for each axis.
dimensions (int):
The number of dimensions in the expansion.
cross_truncation (float, Tuple[float, float]):
Use hyperbolic cross truncation scheme to reduce the number of
terms in expansion. If two values are provided, first is low bound
truncation, while the latter upper bound. If only one value, upper
bound is assumed.
graded (bool):
Graded sorting, meaning the indices are always sorted by the index
sum. E.g. ``(2, 2, 2)`` has a sum of 6, and will therefore be
consider larger than both ``(3, 1, 1)`` and ``(1, 1, 3)``.
reverse (bool):
Reversed lexicographical sorting meaning that ``(1, 3)`` is
considered smaller than ``(3, 1)``, instead of the opposite.
Returns:
list:
Order list of indices.
Examples:
>>> numpoly.glexindex(4).tolist()
[[0], [1], [2], [3]]
>>> numpoly.glexindex(2, dimensions=2).tolist()
[[0, 0], [1, 0], [0, 1]]
>>> numpoly.glexindex(start=2, stop=3, dimensions=2).tolist()
[[2, 0], [1, 1], [0, 2]]
>>> numpoly.glexindex([1, 2, 3]).tolist()
[[0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 2]]
>>> numpoly.glexindex([1, 2, 3], cross_truncation=numpy.inf).tolist()
[[0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 1, 1], [0, 0, 2], [0, 1, 2]]
"""
if stop is None:
start, stop = 0, start
start = numpy.array(start, dtype=int).flatten()
stop = numpy.array(stop, dtype=int).flatten()
start, stop, _ = numpy.broadcast_arrays(start, stop, numpy.empty(dimensions))
cross_truncation = cross_truncation*numpy.ones(2)
indices = _glexindex(start, stop, cross_truncation)
if indices.size:
indices = indices[glexsort(indices.T, graded=graded, reverse=reverse)]
return indices
def _glexindex(start, stop, cross_truncation=1.):
"""Backend for the glexindex function."""
# At the beginning the current list of indices just ranges over the
# last dimension.
bound = stop.max()
dimensions = len(start)
start = numpy.clip(start, a_min=0, a_max=None)
dtype = numpy.uint8 if bound < 256 else numpy.uint16
range_ = numpy.arange(bound, dtype=dtype)
indices = range_[:, numpy.newaxis]
for idx in range(dimensions-1):
# Truncate at each step to keep memory usage low
if idx:
indices = indices[cross_truncate(indices, bound-1, cross_truncation[1])]
# Repeats the current set of indices.
# e.g. [0,1,2] -> [0,1,2,0,1,2,...,0,1,2]
indices = numpy.tile(indices, (bound, 1))
# Stretches ranges over the new dimension.
# e.g. [0,1,2] -> [0,0,...,0,1,1,...,1,2,2,...,2]
front = range_.repeat(len(indices)//bound)[:, numpy.newaxis]
# Puts them two together.
indices = numpy.column_stack((front, indices))
# Complete the truncation scheme
if dimensions == 1:
indices = indices[(indices >= start) & (indices < bound)]
else:
lower = cross_truncate(indices, start-1, cross_truncation[0])
upper = cross_truncate(indices, stop-1, cross_truncation[1])
indices = indices[lower ^ upper]
return numpy.array(indices, dtype=int).reshape(-1, dimensions)
def cross_truncate(indices, bound, norm):
r"""
Truncate of indices using L_p norm.
.. math:
L_p(x) = \sum_i |x_i/b_i|^p ^{1/p} \leq 1
where :math:`b_i` are bounds that each :math:`x_i` should follow.
Args:
indices (Sequence[int]):
Indices to be truncated.
bound (int, Sequence[int]):
The bound function for witch the indices can not be larger than.
norm (float, Sequence[float]):
The `p` in the `L_p`-norm. Support includes both `L_0` and `L_inf`.
Returns:
Boolean indices to ``indices`` with True for each index where the
truncation criteria holds.
Examples:
>>> indices = numpy.array(numpy.mgrid[:10, :10]).reshape(2, -1).T
>>> indices[cross_truncate(indices, 2, norm=0)].T
array([[0, 0, 0, 1, 2],
[0, 1, 2, 0, 0]])
>>> indices[cross_truncate(indices, 2, norm=1)].T
array([[0, 0, 0, 1, 1, 2],
[0, 1, 2, 0, 1, 0]])
>>> indices[cross_truncate(indices, [0, 1], norm=1)].T
array([[0, 0],
[0, 1]])
"""
assert norm >= 0, "negative L_p norm not allowed"
bound = numpy.asfarray(bound).flatten()*numpy.ones(indices.shape[1])
if numpy.any(bound < 0):
return numpy.zeros((len(indices),), dtype=bool)
if numpy.any(bound == 0):
out = numpy.all(indices[:, bound == 0] == 0, axis=-1)
if numpy.any(bound):
out &= cross_truncate(indices[:, bound != 0], bound[bound != 0], norm=norm)
return out
if norm == 0:
out = numpy.sum(indices > 0, axis=-1) <= 1
out[numpy.any(indices > bound, axis=-1)] = False
elif norm == numpy.inf:
out = numpy.max(indices/bound, axis=-1) <= 1
else:
out = numpy.sum((indices/bound)**norm, axis=-1)**(1./norm) <= 1
assert numpy.all(out[numpy.all(indices == 0, axis=-1)])
return out
def glexsort(
keys: numpy.typing.ArrayLike,
graded: bool = False,
reverse: bool = False,
) -> numpy.ndarray:
"""
Sort keys using graded lexicographical ordering.
Same as ``numpy.lexsort``, but also support graded and reverse
lexicographical ordering.
Args:
keys:
Values to sort.
graded:
Graded sorting, meaning the indices are always sorted by the index
sum. E.g. ``(2, 2, 2)`` has a sum of 6, and will therefore be
consider larger than both ``(3, 1, 1)`` and ``(1, 1, 3)``.
reverse:
Reverse lexicographical sorting meaning that ``(1, 3)`` is
considered smaller than ``(3, 1)``, instead of the opposite.
Returns:
Array of indices that sort the keys along the specified axis.
Examples:
>>> indices = numpy.array([[0, 0, 0, 1, 2, 1],
... [1, 2, 0, 0, 0, 1]])
>>> indices[:, numpy.lexsort(indices)]
array([[0, 1, 2, 0, 1, 0],
[0, 0, 0, 1, 1, 2]])
>>> indices[:, numpoly.glexsort(indices)]
array([[0, 1, 2, 0, 1, 0],
[0, 0, 0, 1, 1, 2]])
>>> indices[:, numpoly.glexsort(indices, reverse=True)]
array([[0, 0, 0, 1, 1, 2],
[0, 1, 2, 0, 1, 0]])
>>> indices[:, numpoly.glexsort(indices, graded=True)]
array([[0, 1, 0, 2, 1, 0],
[0, 0, 1, 0, 1, 2]])
>>> indices[:, numpoly.glexsort(indices, graded=True, reverse=True)]
array([[0, 0, 1, 0, 1, 2],
[0, 1, 0, 2, 1, 0]])
>>> indices = numpy.array([4, 5, 6, 3, 2, 1])
>>> indices[numpoly.glexsort(indices)]
array([1, 2, 3, 4, 5, 6])
"""
keys_ = numpy.atleast_2d(keys)
if reverse:
keys_ = keys_[::-1]
indices = numpy.array(numpy.lexsort(keys_))
if graded:
indices = indices[numpy.argsort(
numpy.sum(keys_[:, indices], axis=0))].T
return indices
Functions
def glexindex(start, stop=None, dimensions=1, cross_truncation=1.0, graded=False, reverse=False)
-
Generate graded lexicographical multi-indices for the monomial exponents.
Args
start (Union[int, numpy.ndarray]): The lower order of the indices. If array of int, counts as lower bound for each axis. stop (Union[int, numpy.ndarray, None]): The maximum shape included. If omitted: stop <- start; start <- 0 If int is provided, set as largest total order. If array of int, set as upper bound for each axis. dimensions (int): The number of dimensions in the expansion. cross_truncation (float, Tuple[float, float]): Use hyperbolic cross truncation scheme to reduce the number of terms in expansion. If two values are provided, first is low bound truncation, while the latter upper bound. If only one value, upper bound is assumed. graded (bool): Graded sorting, meaning the indices are always sorted by the index sum. E.g.
(2, 2, 2)
has a sum of 6, and will therefore be consider larger than both(3, 1, 1)
and(1, 1, 3)
. reverse (bool): Reversed lexicographical sorting meaning that(1, 3)
is considered smaller than(3, 1)
, instead of the opposite.Returns
list: Order list of indices.
Examples
>>> numpoly.glexindex(4).tolist() [[0], [1], [2], [3]] >>> numpoly.glexindex(2, dimensions=2).tolist() [[0, 0], [1, 0], [0, 1]] >>> numpoly.glexindex(start=2, stop=3, dimensions=2).tolist() [[2, 0], [1, 1], [0, 2]] >>> numpoly.glexindex([1, 2, 3]).tolist() [[0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 2]] >>> numpoly.glexindex([1, 2, 3], cross_truncation=numpy.inf).tolist() [[0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 1, 1], [0, 0, 2], [0, 1, 2]]
Expand source code
def glexindex(start, stop=None, dimensions=1, cross_truncation=1., graded=False, reverse=False): """ Generate graded lexicographical multi-indices for the monomial exponents. Args: start (Union[int, numpy.ndarray]): The lower order of the indices. If array of int, counts as lower bound for each axis. stop (Union[int, numpy.ndarray, None]): The maximum shape included. If omitted: stop <- start; start <- 0 If int is provided, set as largest total order. If array of int, set as upper bound for each axis. dimensions (int): The number of dimensions in the expansion. cross_truncation (float, Tuple[float, float]): Use hyperbolic cross truncation scheme to reduce the number of terms in expansion. If two values are provided, first is low bound truncation, while the latter upper bound. If only one value, upper bound is assumed. graded (bool): Graded sorting, meaning the indices are always sorted by the index sum. E.g. ``(2, 2, 2)`` has a sum of 6, and will therefore be consider larger than both ``(3, 1, 1)`` and ``(1, 1, 3)``. reverse (bool): Reversed lexicographical sorting meaning that ``(1, 3)`` is considered smaller than ``(3, 1)``, instead of the opposite. Returns: list: Order list of indices. Examples: >>> numpoly.glexindex(4).tolist() [[0], [1], [2], [3]] >>> numpoly.glexindex(2, dimensions=2).tolist() [[0, 0], [1, 0], [0, 1]] >>> numpoly.glexindex(start=2, stop=3, dimensions=2).tolist() [[2, 0], [1, 1], [0, 2]] >>> numpoly.glexindex([1, 2, 3]).tolist() [[0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 2]] >>> numpoly.glexindex([1, 2, 3], cross_truncation=numpy.inf).tolist() [[0, 0, 0], [0, 1, 0], [0, 0, 1], [0, 1, 1], [0, 0, 2], [0, 1, 2]] """ if stop is None: start, stop = 0, start start = numpy.array(start, dtype=int).flatten() stop = numpy.array(stop, dtype=int).flatten() start, stop, _ = numpy.broadcast_arrays(start, stop, numpy.empty(dimensions)) cross_truncation = cross_truncation*numpy.ones(2) indices = _glexindex(start, stop, cross_truncation) if indices.size: indices = indices[glexsort(indices.T, graded=graded, reverse=reverse)] return indices
def cross_truncate(indices, bound, norm)
-
Truncate of indices using L_p norm. .. math: L_p(x) = \sum_i |x_i/b_i|^p ^{1/p} \leq 1 where :math:
b_i
are bounds that each :math:x_i
should follow.Args
indices (Sequence[int]): Indices to be truncated. bound (int, Sequence[int]): The bound function for witch the indices can not be larger than. norm (float, Sequence[float]): The
p
in theL_p
-norm. Support includes bothL_0
andL_inf
.Returns
Boolean indices to
indices
with True for each index where the truncation criteria holds.Examples
>>> indices = numpy.array(numpy.mgrid[:10, :10]).reshape(2, -1).T >>> indices[cross_truncate(indices, 2, norm=0)].T array([[0, 0, 0, 1, 2], [0, 1, 2, 0, 0]]) >>> indices[cross_truncate(indices, 2, norm=1)].T array([[0, 0, 0, 1, 1, 2], [0, 1, 2, 0, 1, 0]]) >>> indices[cross_truncate(indices, [0, 1], norm=1)].T array([[0, 0], [0, 1]])
Expand source code
def cross_truncate(indices, bound, norm): r""" Truncate of indices using L_p norm. .. math: L_p(x) = \sum_i |x_i/b_i|^p ^{1/p} \leq 1 where :math:`b_i` are bounds that each :math:`x_i` should follow. Args: indices (Sequence[int]): Indices to be truncated. bound (int, Sequence[int]): The bound function for witch the indices can not be larger than. norm (float, Sequence[float]): The `p` in the `L_p`-norm. Support includes both `L_0` and `L_inf`. Returns: Boolean indices to ``indices`` with True for each index where the truncation criteria holds. Examples: >>> indices = numpy.array(numpy.mgrid[:10, :10]).reshape(2, -1).T >>> indices[cross_truncate(indices, 2, norm=0)].T array([[0, 0, 0, 1, 2], [0, 1, 2, 0, 0]]) >>> indices[cross_truncate(indices, 2, norm=1)].T array([[0, 0, 0, 1, 1, 2], [0, 1, 2, 0, 1, 0]]) >>> indices[cross_truncate(indices, [0, 1], norm=1)].T array([[0, 0], [0, 1]]) """ assert norm >= 0, "negative L_p norm not allowed" bound = numpy.asfarray(bound).flatten()*numpy.ones(indices.shape[1]) if numpy.any(bound < 0): return numpy.zeros((len(indices),), dtype=bool) if numpy.any(bound == 0): out = numpy.all(indices[:, bound == 0] == 0, axis=-1) if numpy.any(bound): out &= cross_truncate(indices[:, bound != 0], bound[bound != 0], norm=norm) return out if norm == 0: out = numpy.sum(indices > 0, axis=-1) <= 1 out[numpy.any(indices > bound, axis=-1)] = False elif norm == numpy.inf: out = numpy.max(indices/bound, axis=-1) <= 1 else: out = numpy.sum((indices/bound)**norm, axis=-1)**(1./norm) <= 1 assert numpy.all(out[numpy.all(indices == 0, axis=-1)]) return out
def glexsort(keys: Union[numpy._array_like._SupportsArray[numpy.dtype], numpy._nested_sequence._NestedSequence[numpy._array_like._SupportsArray[numpy.dtype]], bool, int, float, complex, str, bytes, numpy._nested_sequence._NestedSequence[Union[bool, int, float, complex, str, bytes]]], graded: bool = False, reverse: bool = False) ‑> numpy.ndarray
-
Sort keys using graded lexicographical ordering. Same as
numpy.lexsort
, but also support graded and reverse lexicographical ordering.Args
keys: Values to sort. graded: Graded sorting, meaning the indices are always sorted by the index sum. E.g.
(2, 2, 2)
has a sum of 6, and will therefore be consider larger than both(3, 1, 1)
and(1, 1, 3)
. reverse: Reverse lexicographical sorting meaning that(1, 3)
is considered smaller than(3, 1)
, instead of the opposite.Returns
Array of indices that sort the keys along the specified axis.
Examples
>>> indices = numpy.array([[0, 0, 0, 1, 2, 1], ... [1, 2, 0, 0, 0, 1]]) >>> indices[:, numpy.lexsort(indices)] array([[0, 1, 2, 0, 1, 0], [0, 0, 0, 1, 1, 2]]) >>> indices[:, numpoly.glexsort(indices)] array([[0, 1, 2, 0, 1, 0], [0, 0, 0, 1, 1, 2]]) >>> indices[:, numpoly.glexsort(indices, reverse=True)] array([[0, 0, 0, 1, 1, 2], [0, 1, 2, 0, 1, 0]]) >>> indices[:, numpoly.glexsort(indices, graded=True)] array([[0, 1, 0, 2, 1, 0], [0, 0, 1, 0, 1, 2]]) >>> indices[:, numpoly.glexsort(indices, graded=True, reverse=True)] array([[0, 0, 1, 0, 1, 2], [0, 1, 0, 2, 1, 0]]) >>> indices = numpy.array([4, 5, 6, 3, 2, 1]) >>> indices[numpoly.glexsort(indices)] array([1, 2, 3, 4, 5, 6])
Expand source code
def glexsort( keys: numpy.typing.ArrayLike, graded: bool = False, reverse: bool = False, ) -> numpy.ndarray: """ Sort keys using graded lexicographical ordering. Same as ``numpy.lexsort``, but also support graded and reverse lexicographical ordering. Args: keys: Values to sort. graded: Graded sorting, meaning the indices are always sorted by the index sum. E.g. ``(2, 2, 2)`` has a sum of 6, and will therefore be consider larger than both ``(3, 1, 1)`` and ``(1, 1, 3)``. reverse: Reverse lexicographical sorting meaning that ``(1, 3)`` is considered smaller than ``(3, 1)``, instead of the opposite. Returns: Array of indices that sort the keys along the specified axis. Examples: >>> indices = numpy.array([[0, 0, 0, 1, 2, 1], ... [1, 2, 0, 0, 0, 1]]) >>> indices[:, numpy.lexsort(indices)] array([[0, 1, 2, 0, 1, 0], [0, 0, 0, 1, 1, 2]]) >>> indices[:, numpoly.glexsort(indices)] array([[0, 1, 2, 0, 1, 0], [0, 0, 0, 1, 1, 2]]) >>> indices[:, numpoly.glexsort(indices, reverse=True)] array([[0, 0, 0, 1, 1, 2], [0, 1, 2, 0, 1, 0]]) >>> indices[:, numpoly.glexsort(indices, graded=True)] array([[0, 1, 0, 2, 1, 0], [0, 0, 1, 0, 1, 2]]) >>> indices[:, numpoly.glexsort(indices, graded=True, reverse=True)] array([[0, 0, 1, 0, 1, 2], [0, 1, 0, 2, 1, 0]]) >>> indices = numpy.array([4, 5, 6, 3, 2, 1]) >>> indices[numpoly.glexsort(indices)] array([1, 2, 3, 4, 5, 6]) """ keys_ = numpy.atleast_2d(keys) if reverse: keys_ = keys_[::-1] indices = numpy.array(numpy.lexsort(keys_)) if graded: indices = indices[numpy.argsort( numpy.sum(keys_[:, indices], axis=0))].T return indices