interval_arithmetic.py 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413
  1. """
  2. Interval Arithmetic for plotting.
  3. This module does not implement interval arithmetic accurately and
  4. hence cannot be used for purposes other than plotting. If you want
  5. to use interval arithmetic, use mpmath's interval arithmetic.
  6. The module implements interval arithmetic using numpy and
  7. python floating points. The rounding up and down is not handled
  8. and hence this is not an accurate implementation of interval
  9. arithmetic.
  10. The module uses numpy for speed which cannot be achieved with mpmath.
  11. """
  12. # Q: Why use numpy? Why not simply use mpmath's interval arithmetic?
  13. # A: mpmath's interval arithmetic simulates a floating point unit
  14. # and hence is slow, while numpy evaluations are orders of magnitude
  15. # faster.
  16. # Q: Why create a separate class for intervals? Why not use SymPy's
  17. # Interval Sets?
  18. # A: The functionalities that will be required for plotting is quite
  19. # different from what Interval Sets implement.
  20. # Q: Why is rounding up and down according to IEEE754 not handled?
  21. # A: It is not possible to do it in both numpy and python. An external
  22. # library has to used, which defeats the whole purpose i.e., speed. Also
  23. # rounding is handled for very few functions in those libraries.
  24. # Q Will my plots be affected?
  25. # A It will not affect most of the plots. The interval arithmetic
  26. # module based suffers the same problems as that of floating point
  27. # arithmetic.
  28. from sympy.core.numbers import int_valued
  29. from sympy.core.logic import fuzzy_and
  30. from sympy.simplify.simplify import nsimplify
  31. from .interval_membership import intervalMembership
  32. class interval:
  33. """ Represents an interval containing floating points as start and
  34. end of the interval
  35. The is_valid variable tracks whether the interval obtained as the
  36. result of the function is in the domain and is continuous.
  37. - True: Represents the interval result of a function is continuous and
  38. in the domain of the function.
  39. - False: The interval argument of the function was not in the domain of
  40. the function, hence the is_valid of the result interval is False
  41. - None: The function was not continuous over the interval or
  42. the function's argument interval is partly in the domain of the
  43. function
  44. A comparison between an interval and a real number, or a
  45. comparison between two intervals may return ``intervalMembership``
  46. of two 3-valued logic values.
  47. """
  48. def __init__(self, *args, is_valid=True, **kwargs):
  49. self.is_valid = is_valid
  50. if len(args) == 1:
  51. if isinstance(args[0], interval):
  52. self.start, self.end = args[0].start, args[0].end
  53. else:
  54. self.start = float(args[0])
  55. self.end = float(args[0])
  56. elif len(args) == 2:
  57. if args[0] < args[1]:
  58. self.start = float(args[0])
  59. self.end = float(args[1])
  60. else:
  61. self.start = float(args[1])
  62. self.end = float(args[0])
  63. else:
  64. raise ValueError("interval takes a maximum of two float values "
  65. "as arguments")
  66. @property
  67. def mid(self):
  68. return (self.start + self.end) / 2.0
  69. @property
  70. def width(self):
  71. return self.end - self.start
  72. def __repr__(self):
  73. return "interval(%f, %f)" % (self.start, self.end)
  74. def __str__(self):
  75. return "[%f, %f]" % (self.start, self.end)
  76. def __lt__(self, other):
  77. if isinstance(other, (int, float)):
  78. if self.end < other:
  79. return intervalMembership(True, self.is_valid)
  80. elif self.start > other:
  81. return intervalMembership(False, self.is_valid)
  82. else:
  83. return intervalMembership(None, self.is_valid)
  84. elif isinstance(other, interval):
  85. valid = fuzzy_and([self.is_valid, other.is_valid])
  86. if self.end < other. start:
  87. return intervalMembership(True, valid)
  88. if self.start > other.end:
  89. return intervalMembership(False, valid)
  90. return intervalMembership(None, valid)
  91. else:
  92. return NotImplemented
  93. def __gt__(self, other):
  94. if isinstance(other, (int, float)):
  95. if self.start > other:
  96. return intervalMembership(True, self.is_valid)
  97. elif self.end < other:
  98. return intervalMembership(False, self.is_valid)
  99. else:
  100. return intervalMembership(None, self.is_valid)
  101. elif isinstance(other, interval):
  102. return other.__lt__(self)
  103. else:
  104. return NotImplemented
  105. def __eq__(self, other):
  106. if isinstance(other, (int, float)):
  107. if self.start == other and self.end == other:
  108. return intervalMembership(True, self.is_valid)
  109. if other in self:
  110. return intervalMembership(None, self.is_valid)
  111. else:
  112. return intervalMembership(False, self.is_valid)
  113. if isinstance(other, interval):
  114. valid = fuzzy_and([self.is_valid, other.is_valid])
  115. if self.start == other.start and self.end == other.end:
  116. return intervalMembership(True, valid)
  117. elif self.__lt__(other)[0] is not None:
  118. return intervalMembership(False, valid)
  119. else:
  120. return intervalMembership(None, valid)
  121. else:
  122. return NotImplemented
  123. def __ne__(self, other):
  124. if isinstance(other, (int, float)):
  125. if self.start == other and self.end == other:
  126. return intervalMembership(False, self.is_valid)
  127. if other in self:
  128. return intervalMembership(None, self.is_valid)
  129. else:
  130. return intervalMembership(True, self.is_valid)
  131. if isinstance(other, interval):
  132. valid = fuzzy_and([self.is_valid, other.is_valid])
  133. if self.start == other.start and self.end == other.end:
  134. return intervalMembership(False, valid)
  135. if not self.__lt__(other)[0] is None:
  136. return intervalMembership(True, valid)
  137. return intervalMembership(None, valid)
  138. else:
  139. return NotImplemented
  140. def __le__(self, other):
  141. if isinstance(other, (int, float)):
  142. if self.end <= other:
  143. return intervalMembership(True, self.is_valid)
  144. if self.start > other:
  145. return intervalMembership(False, self.is_valid)
  146. else:
  147. return intervalMembership(None, self.is_valid)
  148. if isinstance(other, interval):
  149. valid = fuzzy_and([self.is_valid, other.is_valid])
  150. if self.end <= other.start:
  151. return intervalMembership(True, valid)
  152. if self.start > other.end:
  153. return intervalMembership(False, valid)
  154. return intervalMembership(None, valid)
  155. else:
  156. return NotImplemented
  157. def __ge__(self, other):
  158. if isinstance(other, (int, float)):
  159. if self.start >= other:
  160. return intervalMembership(True, self.is_valid)
  161. elif self.end < other:
  162. return intervalMembership(False, self.is_valid)
  163. else:
  164. return intervalMembership(None, self.is_valid)
  165. elif isinstance(other, interval):
  166. return other.__le__(self)
  167. def __add__(self, other):
  168. if isinstance(other, (int, float)):
  169. if self.is_valid:
  170. return interval(self.start + other, self.end + other)
  171. else:
  172. start = self.start + other
  173. end = self.end + other
  174. return interval(start, end, is_valid=self.is_valid)
  175. elif isinstance(other, interval):
  176. start = self.start + other.start
  177. end = self.end + other.end
  178. valid = fuzzy_and([self.is_valid, other.is_valid])
  179. return interval(start, end, is_valid=valid)
  180. else:
  181. return NotImplemented
  182. __radd__ = __add__
  183. def __sub__(self, other):
  184. if isinstance(other, (int, float)):
  185. start = self.start - other
  186. end = self.end - other
  187. return interval(start, end, is_valid=self.is_valid)
  188. elif isinstance(other, interval):
  189. start = self.start - other.end
  190. end = self.end - other.start
  191. valid = fuzzy_and([self.is_valid, other.is_valid])
  192. return interval(start, end, is_valid=valid)
  193. else:
  194. return NotImplemented
  195. def __rsub__(self, other):
  196. if isinstance(other, (int, float)):
  197. start = other - self.end
  198. end = other - self.start
  199. return interval(start, end, is_valid=self.is_valid)
  200. elif isinstance(other, interval):
  201. return other.__sub__(self)
  202. else:
  203. return NotImplemented
  204. def __neg__(self):
  205. if self.is_valid:
  206. return interval(-self.end, -self.start)
  207. else:
  208. return interval(-self.end, -self.start, is_valid=self.is_valid)
  209. def __mul__(self, other):
  210. if isinstance(other, interval):
  211. if self.is_valid is False or other.is_valid is False:
  212. return interval(-float('inf'), float('inf'), is_valid=False)
  213. elif self.is_valid is None or other.is_valid is None:
  214. return interval(-float('inf'), float('inf'), is_valid=None)
  215. else:
  216. inters = []
  217. inters.append(self.start * other.start)
  218. inters.append(self.end * other.start)
  219. inters.append(self.start * other.end)
  220. inters.append(self.end * other.end)
  221. start = min(inters)
  222. end = max(inters)
  223. return interval(start, end)
  224. elif isinstance(other, (int, float)):
  225. return interval(self.start*other, self.end*other, is_valid=self.is_valid)
  226. else:
  227. return NotImplemented
  228. __rmul__ = __mul__
  229. def __contains__(self, other):
  230. if isinstance(other, (int, float)):
  231. return self.start <= other and self.end >= other
  232. else:
  233. return self.start <= other.start and other.end <= self.end
  234. def __rtruediv__(self, other):
  235. if isinstance(other, (int, float)):
  236. other = interval(other)
  237. return other.__truediv__(self)
  238. elif isinstance(other, interval):
  239. return other.__truediv__(self)
  240. else:
  241. return NotImplemented
  242. def __truediv__(self, other):
  243. # Both None and False are handled
  244. if not self.is_valid:
  245. # Don't divide as the value is not valid
  246. return interval(-float('inf'), float('inf'), is_valid=self.is_valid)
  247. if isinstance(other, (int, float)):
  248. if other == 0:
  249. # Divide by zero encountered. valid nowhere
  250. return interval(-float('inf'), float('inf'), is_valid=False)
  251. else:
  252. return interval(self.start / other, self.end / other)
  253. elif isinstance(other, interval):
  254. if other.is_valid is False or self.is_valid is False:
  255. return interval(-float('inf'), float('inf'), is_valid=False)
  256. elif other.is_valid is None or self.is_valid is None:
  257. return interval(-float('inf'), float('inf'), is_valid=None)
  258. else:
  259. # denominator contains both signs, i.e. being divided by zero
  260. # return the whole real line with is_valid = None
  261. if 0 in other:
  262. return interval(-float('inf'), float('inf'), is_valid=None)
  263. # denominator negative
  264. this = self
  265. if other.end < 0:
  266. this = -this
  267. other = -other
  268. # denominator positive
  269. inters = []
  270. inters.append(this.start / other.start)
  271. inters.append(this.end / other.start)
  272. inters.append(this.start / other.end)
  273. inters.append(this.end / other.end)
  274. start = max(inters)
  275. end = min(inters)
  276. return interval(start, end)
  277. else:
  278. return NotImplemented
  279. def __pow__(self, other):
  280. # Implements only power to an integer.
  281. from .lib_interval import exp, log
  282. if not self.is_valid:
  283. return self
  284. if isinstance(other, interval):
  285. return exp(other * log(self))
  286. elif isinstance(other, (float, int)):
  287. if other < 0:
  288. return 1 / self.__pow__(abs(other))
  289. else:
  290. if int_valued(other):
  291. return _pow_int(self, other)
  292. else:
  293. return _pow_float(self, other)
  294. else:
  295. return NotImplemented
  296. def __rpow__(self, other):
  297. if isinstance(other, (float, int)):
  298. if not self.is_valid:
  299. #Don't do anything
  300. return self
  301. elif other < 0:
  302. if self.width > 0:
  303. return interval(-float('inf'), float('inf'), is_valid=False)
  304. else:
  305. power_rational = nsimplify(self.start)
  306. num, denom = power_rational.as_numer_denom()
  307. if denom % 2 == 0:
  308. return interval(-float('inf'), float('inf'),
  309. is_valid=False)
  310. else:
  311. start = -abs(other)**self.start
  312. end = start
  313. return interval(start, end)
  314. else:
  315. return interval(other**self.start, other**self.end)
  316. elif isinstance(other, interval):
  317. return other.__pow__(self)
  318. else:
  319. return NotImplemented
  320. def __hash__(self):
  321. return hash((self.is_valid, self.start, self.end))
  322. def _pow_float(inter, power):
  323. """Evaluates an interval raised to a floating point."""
  324. power_rational = nsimplify(power)
  325. num, denom = power_rational.as_numer_denom()
  326. if num % 2 == 0:
  327. start = abs(inter.start)**power
  328. end = abs(inter.end)**power
  329. if start < 0:
  330. ret = interval(0, max(start, end))
  331. else:
  332. ret = interval(start, end)
  333. return ret
  334. elif denom % 2 == 0:
  335. if inter.end < 0:
  336. return interval(-float('inf'), float('inf'), is_valid=False)
  337. elif inter.start < 0:
  338. return interval(0, inter.end**power, is_valid=None)
  339. else:
  340. return interval(inter.start**power, inter.end**power)
  341. else:
  342. if inter.start < 0:
  343. start = -abs(inter.start)**power
  344. else:
  345. start = inter.start**power
  346. if inter.end < 0:
  347. end = -abs(inter.end)**power
  348. else:
  349. end = inter.end**power
  350. return interval(start, end, is_valid=inter.is_valid)
  351. def _pow_int(inter, power):
  352. """Evaluates an interval raised to an integer power"""
  353. power = int(power)
  354. if power & 1:
  355. return interval(inter.start**power, inter.end**power)
  356. else:
  357. if inter.start < 0 and inter.end > 0:
  358. start = 0
  359. end = max(inter.start**power, inter.end**power)
  360. return interval(start, end)
  361. else:
  362. return interval(inter.start**power, inter.end**power)