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

1"""Convolution of :class:`~chebpy.chebfun.Chebfun` objects. 

2 

3This module implements ``h = f ★ g`` for piecewise Chebfuns. It is kept 

4separate from :mod:`chebpy.chebfun` so the (sizeable) convolution algorithm 

5does not bloat the main class; :meth:`chebpy.chebfun.Chebfun.conv` is a thin 

6wrapper around :func:`convolve`. 

7 

8Two strategies are used: 

9 

10* When both inputs are single-piece :class:`~chebpy.bndfun.Bndfun` funs of 

11 equal width, the fast Hale-Townsend Legendre convolution is used 

12 (:func:`_equal_width_pair`). 

13* Otherwise each output sub-interval is built adaptively using 

14 Gauss-Legendre quadrature (:func:`_piecewise`). 

15 

16The algorithm is based on: 

17 N. Hale and A. Townsend, "An algorithm for the convolution of Legendre 

18 series", SIAM J. Sci. Comput., 36(3), A1207-A1220, 2014. 

19""" 

20 

21from __future__ import annotations 

22 

23from collections.abc import Callable 

24from typing import TYPE_CHECKING, Any 

25 

26import numpy as np 

27 

28from .algorithms import _conv_legendre, cheb2leg, leg2cheb 

29from .bndfun import Bndfun 

30from .chebtech import Chebtech 

31from .fun import Fun 

32from .trigtech import Trigtech 

33from .utilities import Interval 

34 

35if TYPE_CHECKING: 

36 from .chebfun import Chebfun 

37 

38 

39def convolve(f: Chebfun, g: Chebfun) -> Chebfun: 

40 """Return the convolution ``h = f ★ g`` as a piecewise Chebfun. 

41 

42 See :meth:`chebpy.chebfun.Chebfun.conv` for the full description; this is 

43 the implementation behind that method. 

44 """ 

45 if f.isempty or g.isempty: 

46 return f.__class__.initempty() 

47 

48 _reject_unsupported(f, g) 

49 _reject_nonzero_tails(f, g) 

50 

51 # Fast path: both single-piece with equal-width finite domains. 

52 if _use_equal_width_fast_path(f, g): 

53 return _equal_width_pair(f, f.funs[0], g.funs[0]) 

54 

55 # General piecewise convolution. 

56 return _piecewise(f, g) 

57 

58 

59def _reject_unsupported(f: Chebfun, g: Chebfun) -> None: 

60 """Reject convolution operands the algorithms cannot handle. 

61 

62 Both the Hale-Townsend Legendre algorithm and the Gauss-Legendre fallback 

63 assume Chebyshev coefficients on an affine map, so :class:`Trigtech`-backed 

64 pieces (Fourier coefficients; periodic convolution is a distinct 

65 ``circconv`` operation) and :class:`~chebpy.singfun.Singfun` pieces 

66 (non-affine clustering map) are refused with a clear error. 

67 

68 Raises: 

69 NotImplementedError: If either operand contains a Trigtech- or 

70 Singfun-backed piece. 

71 """ 

72 from .singfun import Singfun 

73 

74 pieces = (*f.funs, *g.funs) 

75 if any(isinstance(piece.onefun, Trigtech) for piece in pieces): 

76 raise NotImplementedError( 

77 "conv() is not supported for trigfun (Trigtech-backed) inputs. " 

78 "Aperiodic convolution and periodic (circular) convolution are " 

79 "distinct operations; a dedicated circconv() for trigfuns is " 

80 "not yet implemented." 

81 ) 

82 if any(isinstance(piece, Singfun) for piece in pieces): 

83 raise NotImplementedError( 

84 "conv() is not supported for Chebfuns containing Singfun pieces " 

85 "(functions with endpoint singularities represented by a non-affine " 

86 "clustering map)." 

87 ) 

88 

89 

90def _reject_nonzero_tails(f: Chebfun, g: Chebfun) -> None: 

91 """Reject convolution when a :class:`CompactFun` piece has a non-zero tail. 

92 

93 Convolution of a function with a non-zero asymptotic limit diverges on an 

94 unbounded interval, so refuse early with a clear error pointing the user at 

95 the algebraic-closure escape hatch. 

96 

97 Raises: 

98 DivergentIntegralError: If any CompactFun piece of either operand has a 

99 non-zero ``tail_left`` or ``tail_right``. 

100 """ 

101 from .compactfun import CompactFun 

102 from .exceptions import DivergentIntegralError 

103 

104 for label, h in (("self", f), ("other", g)): 

105 for piece in h.funs: 

106 if isinstance(piece, CompactFun) and (piece.tail_left != 0.0 or piece.tail_right != 0.0): 

107 raise DivergentIntegralError( # noqa: TRY003 

108 f"Convolution requires both operands to decay to zero at " 

109 f"±inf; got tail_left={piece.tail_left}, " 

110 f"tail_right={piece.tail_right} for {label}. Consider " 

111 f"subtracting a matched sigmoid first so the residual " 

112 f"has zero tails, then convolving the residual." 

113 ) 

114 

115 

116def _use_equal_width_fast_path(f: Chebfun, g: Chebfun) -> bool: 

117 """Return True if the fast equal-width Legendre path applies. 

118 

119 The fast path requires both operands to be single finite :class:`Bndfun` 

120 pieces of equal width. :class:`CompactFun` pieces (with possibly infinite 

121 logical support) always take the general piecewise path so the output is 

122 wrapped correctly. 

123 """ 

124 from .compactfun import CompactFun 

125 

126 if any(isinstance(piece, CompactFun) for piece in (*f.funs, *g.funs)): 

127 return False 

128 if f.funs.size != 1 or g.funs.size != 1: 

129 return False 

130 f_fun, g_fun = f.funs[0], g.funs[0] 

131 f_w = float(f_fun.support[1]) - float(f_fun.support[0]) 

132 g_w = float(g_fun.support[1]) - float(g_fun.support[0]) 

133 return bool(np.isclose(f_w, g_w)) 

134 

135 

136def _equal_width_pair(f: Chebfun, f_fun: Any, g_fun: Any) -> Chebfun: 

137 """Convolve two single Bndfuns of equal width using the fast algorithm. 

138 

139 Uses the Hale-Townsend Legendre convolution. The two funs may be on 

140 different intervals as long as they have the same width. 

141 """ 

142 a = float(f_fun.support[0]) 

143 b = float(f_fun.support[1]) 

144 c = float(g_fun.support[0]) 

145 d = float(g_fun.support[1]) 

146 

147 h = (b - a) / 2.0 # half-width (same for both funs) 

148 

149 leg_f = cheb2leg(f_fun.coeffs) 

150 leg_g = cheb2leg(g_fun.coeffs) 

151 

152 gamma_left, gamma_right = _conv_legendre(leg_f, leg_g) 

153 

154 gamma_left = h * gamma_left 

155 gamma_right = h * gamma_right 

156 

157 cheb_left = leg2cheb(gamma_left) 

158 cheb_right = leg2cheb(gamma_right) 

159 

160 mid = (a + b + c + d) / 2.0 

161 left_interval = Interval(a + c, mid) 

162 right_interval = Interval(mid, b + d) 

163 

164 left_fun = Bndfun(Chebtech(cheb_left), left_interval) 

165 right_fun = Bndfun(Chebtech(cheb_right), right_interval) 

166 

167 return f.__class__([left_fun, right_fun]) 

168 

169 

170def _piecewise(f: Chebfun, g: Chebfun) -> Chebfun: 

171 """General piecewise convolution via Gauss-Legendre quadrature. 

172 

173 The breakpoints of the result are the sorted, unique pairwise sums of the 

174 breakpoints of ``f`` and ``g``. On each sub-interval the convolution 

175 integral is smooth, so we construct it adaptively. When either input 

176 contains :class:`CompactFun` pieces, the corresponding ``±inf`` breakpoints 

177 are replaced with the numerical-support bounds for the purposes of 

178 integration; the outermost output pieces are then wrapped as 

179 :class:`CompactFun` so the result preserves the unbounded logical support. 

180 """ 

181 f_logical_breaks = np.array(f.breakpoints, dtype=float) 

182 g_logical_breaks = np.array(g.breakpoints, dtype=float) 

183 left_inf = (not np.isfinite(f_logical_breaks[0])) or (not np.isfinite(g_logical_breaks[0])) 

184 right_inf = (not np.isfinite(f_logical_breaks[-1])) or (not np.isfinite(g_logical_breaks[-1])) 

185 

186 f_breaks = _effective_breakpoints(f) 

187 g_breaks = _effective_breakpoints(g) 

188 

189 # Output breakpoints: all pairwise sums, uniquified and coalesced. 

190 out_breaks = np.unique(np.add.outer(f_breaks, g_breaks).ravel()) 

191 hscl = max(abs(out_breaks[0]), abs(out_breaks[-1]), 1.0) 

192 tol = 10.0 * np.finfo(float).eps * hscl 

193 mask = np.concatenate(([True], np.diff(out_breaks) > tol)) 

194 out_breaks = out_breaks[mask] 

195 

196 conv_eval = _make_evaluator(f, g, f_breaks, g_breaks) 

197 return _build_pieces(f, out_breaks, conv_eval, left_inf=left_inf, right_inf=right_inf) 

198 

199 

200def _effective_breakpoints(h: Chebfun) -> np.ndarray: 

201 """Return ``h``'s breakpoints with ±inf replaced by numerical-support bounds.""" 

202 from .compactfun import CompactFun 

203 

204 bps = np.array(h.breakpoints, dtype=float) 

205 if not np.isfinite(bps[0]) and isinstance(h.funs[0], CompactFun): 

206 bps[0] = float(h.funs[0].numerical_support[0]) 

207 if not np.isfinite(bps[-1]) and isinstance(h.funs[-1], CompactFun): 

208 bps[-1] = float(h.funs[-1].numerical_support[1]) 

209 return bps 

210 

211 

212def _make_evaluator( 

213 f: Chebfun, g: Chebfun, f_breaks: np.ndarray, g_breaks: np.ndarray 

214) -> Callable[[np.ndarray], np.ndarray]: 

215 """Build the Gauss-Legendre quadrature evaluator for ``(f ★ g)``. 

216 

217 The integrand is broken at the breakpoints of ``f`` and the shifted 

218 breakpoints of ``g`` so it is polynomial on each sub-interval; the 

219 quadrature order is chosen to integrate that polynomial exactly. 

220 """ 

221 f_a, f_b = float(f_breaks[0]), float(f_breaks[-1]) 

222 g_c, g_d = float(g_breaks[0]), float(g_breaks[-1]) 

223 

224 max_deg = max(fun.size for fun in f.funs) + max(fun.size for fun in g.funs) 

225 n_quad = max(int(np.ceil((max_deg + 1) / 2)), 16) 

226 quad_nodes, quad_weights = np.polynomial.legendre.leggauss(n_quad) 

227 

228 f_bps = [float(bp) for bp in f_breaks] 

229 g_bps = [float(bp) for bp in g_breaks] 

230 

231 def conv_eval(x: np.ndarray) -> np.ndarray: 

232 """Evaluate (f ★ g)(x) via Gauss-Legendre quadrature.""" 

233 x = np.atleast_1d(np.asarray(x, dtype=float)) 

234 result = np.zeros(x.shape) 

235 for idx in range(x.size): 

236 xi = x[idx] 

237 t_lo = max(f_a, xi - g_d) 

238 t_hi = min(f_b, xi - g_c) 

239 if t_hi <= t_lo: 

240 continue 

241 inner = _subinterval_nodes(xi, t_lo, t_hi, f_bps, g_bps) 

242 total = 0.0 

243 for j in range(len(inner) - 1): 

244 a_int, b_int = inner[j], inner[j + 1] 

245 hw = (b_int - a_int) / 2.0 

246 mid = (a_int + b_int) / 2.0 

247 nodes = hw * quad_nodes + mid 

248 wts = hw * quad_weights 

249 total += np.dot(wts, f(nodes) * g(xi - nodes)) 

250 result[idx] = total 

251 return result 

252 

253 return conv_eval 

254 

255 

256def _subinterval_nodes(xi: float, t_lo: float, t_hi: float, f_bps: list[float], g_bps: list[float]) -> list[float]: 

257 """Return the sorted integration sub-interval boundaries in ``(t_lo, t_hi)``. 

258 

259 Breaks at the breakpoints of ``f`` and the shifted breakpoints of ``g`` 

260 that fall strictly inside ``(t_lo, t_hi)``, keeping the integrand polynomial 

261 on each resulting sub-interval. 

262 """ 

263 inner = [t_lo, t_hi] 

264 inner.extend(bp for bp in f_bps if t_lo < bp < t_hi) 

265 inner.extend(xi - bp for bp in g_bps if t_lo < xi - bp < t_hi) 

266 return sorted(set(inner)) 

267 

268 

269def _build_pieces( 

270 f: Chebfun, 

271 out_breaks: np.ndarray, 

272 conv_eval: Callable[[np.ndarray], np.ndarray], 

273 *, 

274 left_inf: bool, 

275 right_inf: bool, 

276) -> Chebfun: 

277 """Build one fun per output sub-interval, wrapping the unbounded ends. 

278 

279 Interior pieces are finite :class:`Bndfun`; the outermost pieces are wrapped 

280 as :class:`CompactFun` when the corresponding logical edge is ``±inf`` so 

281 the result preserves the unbounded logical support. 

282 """ 

283 from .compactfun import CompactFun 

284 

285 n_pieces = len(out_breaks) - 1 

286 funs_list: list[Fun] = [] 

287 for i in range(n_pieces): 

288 a_storage = float(out_breaks[i]) 

289 b_storage = float(out_breaks[i + 1]) 

290 interval = Interval(a_storage, b_storage) 

291 bnd = Bndfun.initfun_adaptive(conv_eval, interval) 

292 wrap_left = i == 0 and left_inf 

293 wrap_right = i == n_pieces - 1 and right_inf 

294 if wrap_left or wrap_right: 

295 a_logical = -np.inf if wrap_left else a_storage 

296 b_logical = np.inf if wrap_right else b_storage 

297 funs_list.append(CompactFun(bnd.onefun, interval, logical_interval=(a_logical, b_logical))) 

298 else: 

299 funs_list.append(bnd) 

300 

301 return f.__class__(funs_list)