11# Based on an adaptive quadrature algorithm by Pedro Gonnet
22
3+ import functools
34from collections import defaultdict
45from fractions import Fraction
56
@@ -141,6 +142,10 @@ def calc_V(x, n):
141142 return np .array (V ).T
142143
143144
145+ def cached_property (f ):
146+ return property (functools .lru_cache (None )(f ))
147+
148+
144149class Coefficients :
145150 def __init__ (self ) -> None :
146151 self .is_set = False
@@ -162,84 +167,52 @@ def __init__(self) -> None:
162167 # Maximum amount of subdivisions
163168 self .ndiv_max = 20
164169
165- def set ( self ):
166- self . is_set = True
170+ @ cached_property
171+ def xi ( self ):
167172 # The Newton polynomials
168173 xi = [- np .cos (np .linspace (0 , np .pi , n )) for n in self .ns ]
169174 # Make `xi` perfectly anti-symmetric, important for splitting the intervals
170- self ._xi = [(row - row [::- 1 ]) / 2 for row in xi ]
171-
172- # Compute the Vandermonde-like matrix and its inverse.
173- self ._V = [calc_V (x , n ) for x , n in zip (xi , self .ns )]
174- self ._V_inv = list (map (scipy .linalg .inv , self ._V ))
175- self ._Vcond = [
176- scipy .linalg .norm (a , 2 ) * scipy .linalg .norm (b , 2 )
177- for a , b in zip (self ._V , self ._V_inv )
178- ]
179-
180- # Compute the shift matrices.
181- self ._T_left , self ._T_right = [
182- self ._V_inv [3 ] @ calc_V ((xi [3 ] + a ) / 2 , self .ns [3 ]) for a in [- 1 , 1 ]
183- ]
175+ return [(row - row [::- 1 ]) / 2 for row in xi ]
184176
185- # set-up the downdate matrix
186- k = np .arange (self .ns [3 ])
187- self ._alpha = np .sqrt ((k + 1 ) ** 2 / (2 * k + 1 ) / (2 * k + 3 ))
188- self ._gamma = np .concatenate (
189- [[0 , 0 ], np .sqrt (k [2 :] ** 2 / (4 * k [2 :] ** 2 - 1 ))]
190- )
191- self ._b_def = calc_bdef (self .ns )
192-
193- @property
194- def xi (self ):
195- if not self .is_set :
196- self .set ()
197- return self ._xi
198-
199- @property
177+ @cached_property
200178 def V (self ):
201- if not self .is_set :
202- self .set ()
203- return self ._V
179+ # Compute the Vandermonde-like matrix and its inverse.
180+ return [calc_V (x , n ) for x , n in zip (self .xi , self .ns )]
204181
205- @property
182+ @cached_property
206183 def V_inv (self ):
207- if not self .is_set :
208- self .set ()
209- return self ._V_inv
184+ # Compute the inverse Vandermonde-like matrix
185+ return list (map (scipy .linalg .inv , self .V ))
210186
211- @property
187+ @cached_property
212188 def Vcond (self ):
213- if not self .is_set :
214- self .set ()
215- return self ._Vcond
189+ return [
190+ scipy .linalg .norm (a , 2 ) * scipy .linalg .norm (b , 2 )
191+ for a , b in zip (self .V , self .V_inv )
192+ ]
216193
217- @property
194+ @cached_property
195+ def k (self ):
196+ return np .arange (self .ns [3 ])
197+
198+ @cached_property
218199 def T_right (self ):
219- if not self .is_set :
220- self .set ()
221- return self ._T_right
200+ return self .V_inv [3 ] @ calc_V ((self .xi [3 ] + 1 ) / 2 , self .ns [3 ])
222201
223- @property
202+ @cached_property
224203 def T_left (self ):
225- if not self .is_set :
226- self .set ()
227- return self ._T_left
204+ return self .V_inv [3 ] @ calc_V ((self .xi [3 ] - 1 ) / 2 , self .ns [3 ])
228205
229- @property
206+ @cached_property
230207 def alpha (self ):
231- if not self .is_set :
232- self .set ()
233- return self ._alpha
208+ return np .sqrt ((self .k + 1 ) ** 2 / (2 * self .k + 1 ) / (2 * self .k + 3 ))
234209
235- @property
210+ @cached_property
236211 def gamma (self ):
237- if not self . is_set :
238- self .set ()
239- return self . _gamma
212+ return np . concatenate (
213+ [[ 0 , 0 ], np . sqrt ( self .k [ 2 :] ** 2 / ( 4 * self . k [ 2 :] ** 2 - 1 ))]
214+ )
240215
241- @property
216+ @cached_property
242217 def b_def (self ):
243- if not self .is_set :
244- self .set ()
245- return self ._b_def
218+ return calc_bdef (self .ns )
0 commit comments