Coverage for src/chebpy/_pointwise.py: 100%

1"""Pointwise and step-function operations on :class:`~chebpy.chebfun.Chebfun`. 

2 

3Kept separate from :mod:`chebpy.chebfun` so the root-splitting machinery behind 

4``abs``, ``sign``, ``ceil``, ``floor`` and pointwise ``maximum``/``minimum`` 

5does not bloat the main class. Each public Chebfun method is a thin wrapper 

6around the corresponding function here. 

7 

8Every function assumes non-empty operands; the wrapping Chebfun methods handle 

9the empty case (via ``@self_empty()``) before delegating. 

10""" 

11 

12from __future__ import annotations 

13 

14from collections.abc import Callable 

15from typing import TYPE_CHECKING, Any 

16 

17import numpy as np 

18 

19from .bndfun import Bndfun 

20from .exceptions import SupportMismatch 

21from .settings import _preferences as prefs 

22 

23if TYPE_CHECKING: 

24 from .chebfun import Chebfun 

25 

26 

27def absolute(f: Chebfun) -> Chebfun: 

28 """Return ``|f|`` by splitting the domain at the roots of ``f``.""" 

29 newdom = f.domain.merge(f.roots()) 

30 funs = [x.absolute() for x in f._break(newdom)] 

31 return f.__class__(funs) 

32 

33 

34def sign(f: Chebfun) -> Chebfun: 

35 """Return the piecewise sign of ``f``. 

36 

37 Splits the domain at the roots of ``f`` and builds a constant piece with 

38 the appropriate sign on each sub-interval; breakpoints at roots are set to 

39 ``0``. 

40 """ 

41 roots = f.roots() 

42 newdom = f.domain.merge(roots) 

43 funs = [] 

44 for fun in f._break(newdom): 

45 mid = fun.support[0] + 0.5 * (fun.support[-1] - fun.support[0]) 

46 s = float(np.sign(float(f(mid)))) 

47 funs.append(Bndfun.initconst(s, fun.interval)) 

48 result = f.__class__(funs) 

49 # Set breakdata: at roots sign is 0, elsewhere use sign of function. 

50 htol = max(1e2 * f.hscale * prefs.eps, prefs.eps) 

51 for bp in result.breakpoints: 

52 if roots.size > 0 and np.any(np.abs(bp - roots) <= htol): 

53 result.breakdata[bp] = 0.0 

54 else: 

55 result.breakdata[bp] = float(np.sign(float(f(bp)))) 

56 return result 

57 

58 

59def ceil(f: Chebfun) -> Chebfun: 

60 """Return the piecewise ceiling of ``f``.""" 

61 return _step_at_integer_crossings(f, np.ceil) 

62 

63 

64def floor(f: Chebfun) -> Chebfun: 

65 """Return the piecewise floor of ``f``.""" 

66 return _step_at_integer_crossings(f, np.floor) 

67 

68 

69def _step_at_integer_crossings(f: Chebfun, rounder: Callable[[Any], Any]) -> Chebfun: 

70 """Build the piecewise-constant ``rounder(f)`` (used by :func:`ceil`/:func:`floor`). 

71 

72 Splits the domain where ``f`` crosses integer values and assigns each piece 

73 the rounded value at its midpoint. 

74 """ 

75 crossings = _integer_crossings(f) 

76 newdom = f.domain.merge(crossings) 

77 funs = [] 

78 for fun in f._break(newdom): 

79 mid = fun.support[0] + 0.5 * (fun.support[-1] - fun.support[0]) 

80 c = float(rounder(float(f(mid)))) 

81 funs.append(Bndfun.initconst(c, fun.interval)) 

82 result = f.__class__(funs) 

83 for bp in result.breakpoints: 

84 result.breakdata[bp] = float(rounder(float(f(bp)))) 

85 return result 

86 

87 

88def _integer_crossings(f: Chebfun) -> np.ndarray: 

89 """Return the x-values where ``f`` crosses an integer value. 

90 

91 Finds the roots of ``f - n`` for each integer ``n`` in the range of ``f``. 

92 """ 

93 all_values = np.concatenate([fun.values() for fun in f]) 

94 lo = int(np.floor(np.min(all_values))) 

95 hi = int(np.ceil(np.max(all_values))) 

96 crossings = [] 

97 for n in range(lo, hi + 1): 

98 shifted = f - n 

99 crossings.extend(shifted.roots().tolist()) 

100 return np.array(crossings) 

101 

102 

103def maximum_minimum(f: Chebfun, other: Chebfun, comparator: Callable[..., bool]) -> Any: 

104 """Return the pointwise maximum or minimum of ``f`` and ``other``. 

105 

106 ``comparator`` selects which operand wins on each sub-interval: 

107 ``operator.ge`` gives the pointwise maximum and ``operator.lt`` the minimum. 

108 The switch points are the roots of ``f - other``. 

109 """ 

110 if f.isempty or other.isempty: 

111 return f.__class__.initempty() 

112 

113 # Find the intersection of domains. 

114 try: 

115 # Try to use union if supports match. 

116 newdom = f.domain.union(other.domain) 

117 except SupportMismatch: 

118 # If supports don't match, find the intersection. 

119 a_min, a_max = f.support 

120 b_min, b_max = other.support 

121 

122 c_min = max(a_min, b_min) 

123 c_max = min(a_max, b_max) 

124 

125 # If there's no intersection, return empty. 

126 if c_min >= c_max: 

127 return f.__class__.initempty() 

128 

129 # Restrict both functions to the intersection and recurse. 

130 f_restricted = f.restrict([c_min, c_max]) 

131 other_restricted = other.restrict([c_min, c_max]) 

132 return maximum_minimum(f_restricted, other_restricted, comparator) 

133 

134 roots = (f - other).roots() 

135 newdom = newdom.merge(roots) 

136 switch = newdom.support.merge(roots) 

137 

138 if switch.size == 0: # pragma: no cover 

139 return f.__class__.initempty() 

140 

141 keys = 0.5 * ((-1) ** np.arange(switch.size - 1) + 1) 

142 if switch.size > 0 and comparator(other(switch[0]), f(switch[0])): 

143 keys = 1 - keys 

144 funs = np.array([]) 

145 for interval, use_self in zip(switch.intervals, keys, strict=False): 

146 subdom = newdom.restrict(interval) 

147 subfun = f.restrict(subdom) if use_self else other.restrict(subdom) 

148 funs = np.append(funs, subfun.funs) 

149 return f.__class__(funs)