|
2 | 2 |
|
3 | 3 | from collections import defaultdict |
4 | 4 | from fractions import Fraction |
| 5 | +from functools import lru_cache |
5 | 6 |
|
6 | 7 | import numpy as np |
7 | 8 | import scipy.linalg |
8 | 9 |
|
9 | | -try: |
10 | | - from functools import cached_property |
11 | | -except ImportError: |
12 | | - from adaptive.utils import cached_property |
13 | | - |
14 | 10 |
|
15 | 11 | def legendre(n): |
16 | 12 | """Return the first n Legendre polynomials. |
@@ -146,72 +142,48 @@ def calc_V(x, n): |
146 | 142 | return np.array(V).T |
147 | 143 |
|
148 | 144 |
|
149 | | -class Coefficients: |
150 | | - def __init__(self) -> None: |
151 | | - # The nodes |
152 | | - self.ns = (5, 9, 17, 33) |
153 | | - |
154 | | - # If the relative difference between two consecutive approximations is |
155 | | - # lower than this value, the error estimate is considered reliable. |
156 | | - # See section 6.2 of Pedro Gonnet's thesis. |
157 | | - self.hint = 0.1 |
158 | | - |
159 | | - # Smallest acceptable relative difference of points in a rule. This was chosen |
160 | | - # such that no artifacts are apparent in plots of (i, log(a_i)), where a_i is |
161 | | - # the sequence of estimates of the integral value of an interval and all its |
162 | | - # ancestors.. |
163 | | - self.eps = np.spacing(1) |
164 | | - self.min_sep = 16 * self.eps |
165 | | - |
166 | | - # Maximum amount of subdivisions |
167 | | - self.ndiv_max = 20 |
168 | | - |
169 | | - @cached_property |
170 | | - def xi(self): |
171 | | - # The Newton polynomials |
172 | | - xi = [-np.cos(np.linspace(0, np.pi, n)) for n in self.ns] |
173 | | - # Make `xi` perfectly anti-symmetric, important for splitting the intervals |
174 | | - return [(row - row[::-1]) / 2 for row in xi] |
175 | | - |
176 | | - @cached_property |
177 | | - def V(self): |
178 | | - # Compute the Vandermonde-like matrix and its inverse. |
179 | | - return [calc_V(x, n) for x, n in zip(self.xi, self.ns)] |
180 | | - |
181 | | - @cached_property |
182 | | - def V_inv(self): |
183 | | - # Compute the inverse Vandermonde-like matrix |
184 | | - return list(map(scipy.linalg.inv, self.V)) |
185 | | - |
186 | | - @cached_property |
187 | | - def Vcond(self): |
188 | | - return [ |
189 | | - scipy.linalg.norm(a, 2) * scipy.linalg.norm(b, 2) |
190 | | - for a, b in zip(self.V, self.V_inv) |
191 | | - ] |
192 | | - |
193 | | - @cached_property |
194 | | - def k(self): |
195 | | - return np.arange(self.ns[3]) |
196 | | - |
197 | | - @cached_property |
198 | | - def T_right(self): |
199 | | - return self.V_inv[3] @ calc_V((self.xi[3] + 1) / 2, self.ns[3]) |
200 | | - |
201 | | - @cached_property |
202 | | - def T_left(self): |
203 | | - return self.V_inv[3] @ calc_V((self.xi[3] - 1) / 2, self.ns[3]) |
204 | | - |
205 | | - @cached_property |
206 | | - def alpha(self): |
207 | | - return np.sqrt((self.k + 1) ** 2 / (2 * self.k + 1) / (2 * self.k + 3)) |
208 | | - |
209 | | - @cached_property |
210 | | - def gamma(self): |
211 | | - return np.concatenate( |
212 | | - [[0, 0], np.sqrt(self.k[2:] ** 2 / (4 * self.k[2:] ** 2 - 1))] |
213 | | - ) |
214 | | - |
215 | | - @cached_property |
216 | | - def b_def(self): |
217 | | - return calc_bdef(self.ns) |
| 145 | +@lru_cache(maxsize=None) |
| 146 | +def _coefficients(): |
| 147 | + eps = np.spacing(1) |
| 148 | + |
| 149 | + # the nodes and Newton polynomials |
| 150 | + ns = (5, 9, 17, 33) |
| 151 | + xi = [-np.cos(np.linspace(0, np.pi, n)) for n in ns] |
| 152 | + |
| 153 | + # Make `xi` perfectly anti-symmetric, important for splitting the intervals |
| 154 | + xi = [(row - row[::-1]) / 2 for row in xi] |
| 155 | + |
| 156 | + # Compute the Vandermonde-like matrix and its inverse. |
| 157 | + V = [calc_V(x, n) for x, n in zip(xi, ns)] |
| 158 | + V_inv = list(map(scipy.linalg.inv, V)) |
| 159 | + Vcond = [ |
| 160 | + scipy.linalg.norm(a, 2) * scipy.linalg.norm(b, 2) for a, b in zip(V, V_inv) |
| 161 | + ] |
| 162 | + |
| 163 | + # Compute the shift matrices. |
| 164 | + T_left, T_right = [V_inv[3] @ calc_V((xi[3] + a) / 2, ns[3]) for a in [-1, 1]] |
| 165 | + |
| 166 | + # If the relative difference between two consecutive approximations is |
| 167 | + # lower than this value, the error estimate is considered reliable. |
| 168 | + # See section 6.2 of Pedro Gonnet's thesis. |
| 169 | + hint = 0.1 |
| 170 | + |
| 171 | + # Smallest acceptable relative difference of points in a rule. This was chosen |
| 172 | + # such that no artifacts are apparent in plots of (i, log(a_i)), where a_i is |
| 173 | + # the sequence of estimates of the integral value of an interval and all its |
| 174 | + # ancestors.. |
| 175 | + min_sep = 16 * eps |
| 176 | + |
| 177 | + ndiv_max = 20 |
| 178 | + |
| 179 | + # set-up the downdate matrix |
| 180 | + k = np.arange(ns[3]) |
| 181 | + alpha = np.sqrt((k + 1) ** 2 / (2 * k + 1) / (2 * k + 3)) |
| 182 | + gamma = np.concatenate([[0, 0], np.sqrt(k[2:] ** 2 / (4 * k[2:] ** 2 - 1))]) |
| 183 | + |
| 184 | + b_def = calc_bdef(ns) |
| 185 | + return locals() |
| 186 | + |
| 187 | + |
| 188 | +def __getattr__(attr): |
| 189 | + return _coefficients()[attr] |
0 commit comments