Coverage for /opt/conda/envs/apienv/lib/python3.10/site-packages/daiquiri/core/transform/sets.py: 92%
249 statements
« prev ^ index » next coverage.py v7.6.5, created at 2024-11-15 02:12 +0000
1"""
2Sets in the mathematical sense
4 Set: any set of elements
5 RealSubset: set of real numbers
6 RealInterval: R^n with limits [a, b] (open, closed or half-open)
7 on each dimension (you can use inf)
8 RealVectorSpace: R^n
9"""
10import sys
11import itertools
12from abc import ABC, abstractmethod, abstractproperty
13import numpy
14from daiquiri.core.transform import utils
16from daiquiri.core.utils import timed
19class DomainError(ValueError):
20 pass
23class NotInSetError(DomainError):
24 pass
27class DimensionError(DomainError):
28 pass
31class Set(ABC):
32 """A set in the mathematical sense"""
34 def __contains__(self, other):
35 """
36 :param num or array-like or Set other: element(s) or set
37 :returns bool:
38 """
39 if isinstance(other, Set):
40 return self.isSubset(other)
41 else:
42 return self.isElement(other)
44 def __eq__(self, other):
45 return other in self and self in other
47 def __ne__(self, other):
48 return not self.__eq__(other)
50 def __le__(self, other):
51 return self in other
53 def __gt__(self, other):
54 return self not in other
56 def __ge__(self, other):
57 return other in self
59 def __lt__(self, other):
60 return other not in self
62 @abstractmethod
63 def isSubset(self, other):
64 """Is other a subset of self (not strict)
66 :param Set other:
67 :returns bool:
68 """
69 pass
71 def isElement(self, elements):
72 """
73 :param num or sequence elements:
74 :returns bool:
75 """
76 return all(self.isElementAsSequence(elements))
78 @abstractmethod
79 def isElementAsSequence(self, elements):
80 """
81 :param num or sequence elements:
82 :returns sequence(bool):
83 """
84 pass
86 @abstractproperty
87 def isFinite(self):
88 """
89 :returns bool:
90 """
91 pass
93 @property
94 def isInfinite(self):
95 """
96 :returns bool:
97 """
98 return not self.isFinite
100 def raiseIfNotIn(self, other):
101 """
102 :param num or array-like or Set other: element(s) or set
103 :returns bool:
104 """
105 if isinstance(other, Set):
106 if other not in self:
107 raise ValueError("Is not a subset")
108 else:
109 self.raiseIfNotIsElement(other)
111 # @timed
112 def raiseIfNotIsElement(self, elements):
113 """
114 :param num or sequence elements:
115 :returns bool:
116 """
117 # TODO: This takes of the order of 10-50ms per call (!)
118 # due to `isElementAsSequence`, disabled for now
119 if not hasattr(sys, "_running_pytest"):
120 return
122 elements = self.as_sequence(elements)
123 inside = self.isElementAsSequence(elements)
124 if not inside.all():
125 outside = elements[~inside]
126 raise NotInSetError(
127 "These elements are not part of the set {}:\n {}".format(
128 repr(self), outside
129 )
130 )
132 def as_sequence(self, elements):
133 """
134 :param num or sequence elements:
135 :returns array:
136 """
137 return numpy.asarray(elements)
140class RealSubset(Set):
141 """Subset of real vector space R^n"""
143 def __init__(self, ndim=None, fillvalue=None):
144 if ndim is None or ndim < 1:
145 raise ValueError("Specify a dimension >= 1")
146 self.__ndim = ndim
147 if fillvalue is None:
148 fillvalue = numpy.nan
149 fillvalue = numpy.atleast_1d(fillvalue)
150 if fillvalue.shape == (1,):
151 fillvalue = numpy.full((ndim,), fillvalue)
152 elif fillvalue.shape != (ndim,):
153 raise ValueError(f"Expected one or {ndim} fill values")
154 self.__fillvalue = fillvalue
156 @property
157 def ndim(self):
158 """Dimension of R^n"""
159 return self.__ndim
161 def as_sequence(self, points):
162 """
163 :param num or sequence points: shape is (), (ndim,) or (npoints, ndim)
164 :returns array: shape is (npoints, ndim)
165 """
166 points = numpy.atleast_2d(points)
167 ndimextra = self.ndim - points.shape[1]
168 if points.ndim != 2 or ndimextra < 0:
169 raise DimensionError(
170 "Coordinate shape must be (npoints, {})".format(self.ndim)
171 )
172 if ndimextra:
173 addnan = self.__fillvalue[numpy.newaxis, -ndimextra:]
174 addnan = numpy.tile(addnan, (points.shape[0], 1))
175 points = numpy.hstack([points, addnan])
176 return points
178 def isSubset(self, other):
179 """Is other a subset of self (not strict)
181 :param RealVectorSpace other:
182 :returns bool:
183 """
184 if not isinstance(other, RealInterval):
185 raise ValueError("Must be an instance of RealInterval")
186 if self.isVectorSpace:
187 return other.ndim <= self.ndim
188 else:
189 return self._isSubset(other)
191 @abstractmethod
192 def _isSubset(self, other):
193 """Is other a subset of self (not strict)
195 :param RealVectorSpace other:
196 :returns bool:
197 """
198 pass
200 def isElementAsSequence(self, points):
201 """
202 :param num or sequence points:
203 :returns array(bool): shape is (npoints,)
204 """
205 return self.isElementAsArray(points).all(axis=1)
207 def isElementAsArray(self, points):
208 """
209 :param num or sequence points:
210 :returns array(bool): shape is (npoints, ndim)
211 """
212 try:
213 points = self.as_sequence(points)
214 except DimensionError:
215 points = numpy.atleast_2d(points)
216 return numpy.full(points.shape, False)
217 return self._isElementAsArray(points)
219 @abstractmethod
220 def _isElementAsArray(self, points):
221 """
222 :param array points: shape is (npoints, ndim)
223 :returns array(bool): shape is (npoints, ndim)
224 """
225 pass
227 @property
228 def isVectorSpace(self):
229 """
230 :returns bool: ndim
231 """
232 return False
235class RealInterval(RealSubset):
236 """Subset of real vector space R^n defined by limits (open, closed, half-open) in each dimensions"""
238 def __init__(self, ndim=None, limits=None, **kwargs):
239 """
240 :param num ndim: dimension of real vector space
241 :param num or array-like limits: limits with shape (ndim, 4) or (ndim, 2)
242 For example [-1, 1, True, False]
243 refers to the interval [-1, 1).
244 """
245 if ndim is None and limits is not None:
246 limits = numpy.atleast_2d(limits)
247 ndim = len(limits)
248 super().__init__(ndim=ndim, **kwargs)
249 self._limits = self._parse_limits(limits)
251 @property
252 def limits(self):
253 """
254 :returns recarray: shape is (ndim, 4)
255 """
256 return self._limits
258 def __repr__(self):
259 if self.isVectorSpace:
260 return "R^{}".format(self.ndim)
261 else:
262 return "R^{}(interval={})".format(self.ndim, self.limits)
264 def _parse_limits(self, values):
265 """Closed by default
267 :param num or array-like values: limits with shape (ndim, 4) or (ndim, 2)
268 :returns recarray: shape is (ndim, 4)
269 """
270 if values is None:
271 values = [[-numpy.inf, numpy.inf]] * self.ndim
272 values = numpy.atleast_2d(values)
273 shape = (self.ndim, 2)
274 if values.shape == shape:
275 cvalues = numpy.full(shape, True)
276 values = numpy.concatenate([values, cvalues], axis=1)
277 shape = (self.ndim, 4)
278 if values.shape != shape:
279 raise DimensionError("Required shape is {}".format(shape))
280 idx = utils.argsort(values[:, :2], axis=1)
281 values[:, :2] = values[:, :2][idx]
282 values[:, 2:] = values[:, 2:][idx]
283 values = numpy.array([tuple(v) for v in values], dtype=self._limit_dtype)
284 return values.view(numpy.recarray)
286 @property
287 def _limit_dtype(self):
288 return [("low", "<f8"), ("high", "<f8"), ("clow", "?"), ("chigh", "?")]
290 def _isSubset(self, other):
291 """Is other a subset of self (not strict)
293 :param RealVectorSpace other:
294 :returns bool:
295 """
296 if not isinstance(other, RealInterval):
297 raise ValueError("Must be an instance of RealInterval")
298 if self.isVectorSpace:
299 return other.ndim <= self.ndim
300 else:
301 return other.border_inlim in self
303 @timed
304 def _isElementAsArray(self, points):
305 """
306 :param array points: shape is (npoints, ndim)
307 :returns array(bool): shape is (npoints, ndim)
308 """
309 # TODO: This takes a long time to execute, why?
310 limits = self.limits
311 b = numpy.full(points.shape, True)
312 if self.isVectorSpace:
313 return b
314 # Lower limit (open or closed)
315 clow = limits.clow
316 olow = ~clow
317 with numpy.errstate(invalid="ignore"):
318 # Ignore nan's
319 b[:, clow] &= points[:, clow] >= limits.low[numpy.newaxis, clow]
320 b[:, olow] &= points[:, olow] > limits.low[numpy.newaxis, olow]
321 # Upper limit (open or closed)
322 chigh = limits.chigh
323 ohigh = ~chigh
324 with numpy.errstate(invalid="ignore"):
325 # Ignore nan's
326 b[:, chigh] &= points[:, chigh] <= limits.high[numpy.newaxis, chigh]
327 b[:, ohigh] &= points[:, ohigh] < limits.high[numpy.newaxis, ohigh]
328 # Nan's: coordinate in subspace
329 b[numpy.isnan(points)] = True
330 return b
332 @property
333 def border(self):
334 """Could be infinite or open
336 :returns array: shape is (npoints, ndim)
337 """
338 return numpy.asarray(
339 list(itertools.product(*zip(self.limits.low, self.limits.high)))
340 )
342 def border_finite(self, maxnum=None):
343 """Could be open limit
345 :returns array: shape is (npoints, ndim)
346 """
347 return numpy.asarray(
348 list(
349 itertools.product(
350 *zip(
351 self._low_finite(maxnum=maxnum),
352 self._high_finite(maxnum=maxnum),
353 )
354 )
355 )
356 )
358 @property
359 def border_inlim(self):
360 """Could be infinite
362 :returns array: shape is (npoints, ndim)
363 """
364 return numpy.asarray(
365 list(itertools.product(*zip(self._low_inlim, self._high_inlim)))
366 )
368 @property
369 def center(self):
370 """Center of the domain
372 :returns array: ndim
373 """
374 cen = (self._high_finite() + self._low_finite()) / 2.0
375 cen[self.fullDimension] = 0
376 return cen
378 @property
379 def length(self):
380 """Length of the interval
382 :returns array: ndim
383 """
384 return self.limits.high - self.limits.low
386 @property
387 def length_inlim(self):
388 """Length of the interval
390 :returns array: ndim
391 """
392 return self._high_inlim - self._low_inlim
394 def random_uniform(self, n, maxnum=None):
395 """
396 :param num n: number of points
397 :param num maxnum:
398 :returns array: shape (n, ndim)
399 """
400 # numpy.random.uniform: [low, high)
402 # Open lower limit when needed
403 low = self._low_finite(maxnum=maxnum)
404 lowopen = ~self.limits.clow
405 low = self._nextafter(low, lowopen, self._high_inlim)
407 # Close higher limit when needed
408 high = self._high_finite(maxnum=maxnum)
409 highclosed = self.limits.chigh & (low != high)
410 high = self._nextafter(high, highclosed, numpy.inf)
411 high = utils.makeFinite(high, maxnum=maxnum)
413 return numpy.random.uniform(low, high, (n, self.ndim))
415 @property
416 def isFinite(self):
417 """
418 :returns bool:
419 """
420 return self.finite_dimensions.all()
422 @property
423 def finite_dimensions(self):
424 """
425 :returns array(bool): ndim
426 """
427 return numpy.isfinite(self.length)
429 @property
430 def fullDimension(self):
431 """
432 :returns array(bool): ndim
433 """
434 return (self.limits.low == -numpy.inf) & (self.limits.high == numpy.inf)
436 @property
437 def isVectorSpace(self):
438 """
439 :returns bool: ndim
440 """
441 return self.fullDimension.all()
443 def _low_finite(self, maxnum=None):
444 """Finite but outside the limits when open
446 :returns array: ndim
447 """
448 return utils.makeFinite(self.limits.low.copy(), maxnum=maxnum)
450 def _high_finite(self, maxnum=None):
451 """Finite but outside the limits when open
453 :returns array: ndim
454 """
455 return utils.makeFinite(self.limits.high.copy(), maxnum=maxnum)
457 def _low_finite_inlim(self, maxnum=None):
458 return utils.makeFinite(self._low_inlim, maxnum=maxnum)
460 def _high_finite_inlim(self, maxnum=None):
461 return utils.makeFinite(self._high_inlim, maxnum=maxnum)
463 @property
464 def _low_inlim(self):
465 """Inside the limits (respect open limit)
467 :returns array: ndim
468 """
469 arr = self.limits.low
470 isopen = ~self.limits.clow
471 towards = self.limits.high
472 return self._nextafter(arr, isopen, towards)
474 @property
475 def _high_inlim(self):
476 """Inside the limits (respect open limit)
478 :returns array: ndim
479 """
480 arr = self.limits.high
481 isopen = ~self.limits.chigh
482 towards = self.limits.low
483 return self._nextafter(arr, isopen, towards)
485 @staticmethod
486 def _nextafter(arr, select, direction):
487 """
488 :param array arr:
489 :param array(bool) select: mask to find next
490 :param num or array direction: direction to find the next dtype value
491 """
492 arr = numpy.tile(arr, (2, 1))
493 try:
494 direction = direction[select]
495 except TypeError:
496 pass
497 arr[1][select] = direction
498 with numpy.errstate(under="ignore", over="ignore"):
499 return numpy.nextafter(*arr)
502class RealIntervalComposite(RealInterval):
503 """Combine several Real intervals"""
505 def __init__(self, intervals, **kwargs):
506 """
507 :param list(RealInterval) intervals:
508 """
509 self._intervals = intervals
510 ndim = sum(d.ndim for d in intervals)
511 super().__init__(ndim=ndim, **kwargs)
513 @property
514 def limits(self):
515 limits = []
516 for interval in self._intervals:
517 limits.extend(interval.limits.tolist())
518 return self._parse_limits(limits)
521class RealVectorSpace(RealInterval):
522 """R^n"""
524 def __init__(self, ndim=None, **kwargs):
525 """
526 :param num ndim: dimension of real vector space
527 """
528 super().__init__(ndim=ndim, **kwargs)