Source code for numdifftools.extrapolation

"""
Created on 28. aug. 2015

@author: pab
"""
from __future__ import division, print_function
import numpy as np
from scipy import linalg
from scipy.ndimage.filters import convolve1d
import warnings

EPS = np.finfo(float).eps
_EPS = EPS
_TINY = np.finfo(float).tiny


[docs]def convolve(sequence, rule, **kwds): """Wrapper around scipy.ndimage.convolve1d that allows complex input.""" if np.iscomplexobj(sequence): return (convolve1d(sequence.real, rule, **kwds) + 1j * convolve1d(sequence.imag, rule, **kwds)) return convolve1d(sequence, rule, **kwds)
[docs]class Dea(object): """ Extrapolate a slowly convergent sequence LIMEXP is the maximum number of elements the epsilon table data can contain. The epsilon table is stored in the first (LIMEXP+2) entries of EPSTAB. LIST OF MAJOR VARIABLES ----------------------- E0,E1,E2,E3 - DOUBLE PRECISION The 4 elements on which the computation of a new element in the epsilon table is based. NRES - INTEGER Number of extrapolation results actually generated by the epsilon algorithm in prior calls to the routine. NEWELM - INTEGER Number of elements to be computed in the new diagonal of the epsilon table. The condensed epsilon table is computed. Only those elements needed for the computation of the next diagonal are preserved. RES - DOUBLE PREISION New element in the new diagonal of the epsilon table. ERROR - DOUBLE PRECISION An estimate of the absolute error of RES. Routine decides whether RESULT=RES or RESULT=SVALUE by comparing ERROR with abserr from the previous call. RES3LA - DOUBLE PREISION Vector of DIMENSION 3 containing at most the last 3 results. """
[docs] def __init__(self, limexp=3): self.limexp = limexp self.abserr = 10. self._n = 0 self._nres = 0
@property def limexp(self): return self._limexp @limexp.setter def limexp(self, limexp): n = 2 * (limexp // 2) + 1 if n < 3: raise ValueError('LIMEXP IS LESS THAN 3') self.epstab = np.zeros(n + 5) self._limexp = n @staticmethod def _compute_error(res3la, nres, res): fact = [6.0, 2.0, 1.0][min(nres - 1, 2)] error = fact * np.abs(res - res3la[:nres]).sum() return error @staticmethod def _shift_table(epstab, n, newelm, old_n): i_0 = old_n % 2 # 1 if ((old_n // 2) * 2 == old_n - 1) else 0 i_n = 2 * newelm + 2 epstab[i_0:i_n:2] = epstab[i_0 + 2:i_n + 2:2] if old_n != n: i_n = old_n - n epstab[:n + 1] = epstab[i_n:i_n + n + 1] return epstab @staticmethod def _update_res3la(res3la, result, nres): if nres > 2: res3la[:2] = res3la[1:] res3la[2] = result else: res3la[nres] = result def _dea(self, epstab, n): nres = self._nres res3la = epstab[-3:] abserr = self.abserr epstab[n + 2] = epstab[n] newelm = n // 2 old_n = n k1 = n for I in range(newelm): e0, e1, e2 = epstab[k1 - 2], epstab[k1 - 1], epstab[k1 + 2] res = e2 delta2, delta3 = e2 - e1, e1 - e0 err2, err3 = abs(delta2), abs(delta3) tol2 = max(abs(e2), abs(e1)) * _EPS tol3 = max(abs(e1), abs(e0)) * _EPS all_converged = err2 <= tol2 and err3 <= tol3 if all_converged: abserr = err2 + err3 result = res break if I == 0: any_converged = err2 <= tol2 or err3 <= tol3 if not any_converged: SS = 1.0 / delta2 - 1.0 / delta3 else: e3 = epstab[k1] delta1 = e1 - e3 err1 = abs(delta1) tol1 = max(abs(e1), abs(e3)) * _EPS any_converged = err1 <= tol1 or err2 <= tol2 or err3 <= tol3 if not any_converged: SS = 1.0 / delta1 + 1.0 / delta2 - 1.0 / delta3 epstab[k1] = e1 if any_converged or abs(SS * e1) <= 1e-04: n = 2 * I if nres == 0: abserr = err2 + err3 result = res else: result = res3la[min(nres - 1, 2)] break res = e1 + 1.0 / SS epstab[k1] = res k1 = k1 - 2 if nres == 0: abserr = err2 + abs(res - e2) + err3 result = res continue error = self._compute_error(res3la, nres, res) if error > 10.0 * abserr: continue abserr = error result = res # else: # pass # error = self._compute_error(res3la, nres, res) # result = res # 50 if n == self.limexp - 1: n = 2 * (self.limexp // 2) - 1 epstab = self._shift_table(epstab, n, newelm, old_n) self._update_res3la(res3la, result, nres) abserr = max(abserr, 10.0 * _EPS * abs(result)) self._nres += 1 return result, abserr, n def __call__(self, s_value): epstab = self.epstab result = s_value n = self._n epstab[n] = s_value if n == 0: abserr = abs(result) elif n == 1: abserr = 6.0 * abs(result - epstab[0]) else: result, abserr, n = self._dea(epstab, n) n += 1 self._n = n self.abserr = abserr return result, abserr
class EpsAlg(object): """ Extrapolate a slowly convergent sequence This implementaion is from [1]_ Reference --------- .. [1] E. J. Weniger (1989) "Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series" Computer Physics Reports Vol. 10, 189 - 371 http://arxiv.org/abs/math/0306302v1 """ def __init__(self, limexp=3): self.limexp = 2 * (limexp // 2) + 1 self.epstab = np.zeros(limexp + 5) self.abserr = 10. self._n = 0 self._nres = 0 if limexp < 3: raise ValueError('LIMEXP IS LESS THAN 3') def __call__(self, s_n): n = self._n epstab = self.epstab epstab[n] = s_n if n == 0: estlim = s_n else: aux2 = 0.0 for i in range(n, 0, -1): aux1 = aux2 aux2 = epstab[i - 1] delta = epstab[i] - aux2 if np.abs(delta) <= 1.0e-60: epstab[i - 1] = 1.0e+60 else: epstab[i - 1] = aux1 + 1.0 / delta estlim = epstab[n % 2] if n > self.limexp - 1: raise ValueError("Eps table to small!") n += 1 self._n = n return estlim def epsalg_demo(): def linfun(i): return np.linspace(0, np.pi / 2., 2**i + 1) dea = EpsAlg(limexp=15) print('NO. PANELS TRAP. APPROX APPROX W/EA abserr') txt = '{0:5d} {1:20.8f} {2:20.8f} {3:20.8f}' for k in np.arange(10): x = linfun(k) val = np.trapz(np.sin(x), x) vale = dea(val) err = np.abs(1.0 - vale) print(txt.format(len(x) - 1, val, vale, err)) def dea_demo(): def linfun(i): return np.linspace(0, np.pi / 2., 2**i + 1) dea = Dea(limexp=6) print('NO. PANELS TRAP. APPROX APPROX W/EA abserr') txt = '{0:5d} {1:20.8f} {2:20.8f} {3:20.8f}' for k in np.arange(12): x = linfun(k) val = np.trapz(np.sin(x), x) vale, err = dea(val) print(txt.format(len(x) - 1, val, vale, err)) def max_abs(a1, a2): return np.maximum(np.abs(a1), np.abs(a2))
[docs]def dea3(v0, v1, v2, symmetric=False): """ Extrapolate a slowly convergent sequence Parameters ---------- v0, v1, v2 : array-like 3 values of a convergent sequence to extrapolate Returns ------- result : array-like extrapolated value abserr : array-like absolute error estimate Description ----------- DEA3 attempts to extrapolate nonlinearly to a better estimate of the sequence's limiting value, thus improving the rate of convergence. The routine is based on the epsilon algorithm of P. Wynn, see [1]_. Example ------- # integrate sin(x) from 0 to pi/2 >>> import numpy as np >>> import numdifftools as nd >>> Ei= np.zeros(3) >>> linfun = lambda i : np.linspace(0, np.pi/2., 2**(i+5)+1) >>> for k in np.arange(3): ... x = linfun(k) ... Ei[k] = np.trapz(np.sin(x),x) >>> [En, err] = nd.dea3(Ei[0], Ei[1], Ei[2]) >>> truErr = Ei-1. >>> (truErr, err, En) (array([ -2.00805680e-04, -5.01999079e-05, -1.25498825e-05]), array([ 0.00020081]), array([ 1.])) See also -------- dea Reference --------- .. [1] C. Brezinski and M. Redivo Zaglia (1991) "Extrapolation Methods. Theory and Practice", North-Holland. .. [2] C. Brezinski (1977) "Acceleration de la convergence en analyse numerique", "Lecture Notes in Math.", vol. 584, Springer-Verlag, New York, 1977. .. [3] E. J. Weniger (1989) "Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series" Computer Physics Reports Vol. 10, 189 - 371 http://arxiv.org/abs/math/0306302v1 """ e0, e1, e2 = np.atleast_1d(v0, v1, v2) with warnings.catch_warnings(): warnings.simplefilter("ignore") # ignore division by zero and overflow delta2, delta1 = e2 - e1, e1 - e0 err2, err1 = np.abs(delta2), np.abs(delta1) tol2, tol1 = max_abs(e2, e1) * _EPS, max_abs(e1, e0) * _EPS delta1[err1 < _TINY] = _TINY delta2[err2 < _TINY] = _TINY # avoid division by zero and overflow ss = 1.0 / delta2 - 1.0 / delta1 + _TINY smalle2 = abs(ss * e1) <= 1.0e-3 converged = (err1 <= tol1) & (err2 <= tol2) | smalle2 result = np.where(converged, e2 * 1.0, e1 + 1.0 / ss) abserr = err1 + err2 + np.where(converged, tol2 * 10, np.abs(result - e2)) if symmetric and len(result) > 1: return result[:-1], abserr[1:] return result, abserr
[docs]class Richardson(object): """ Extrapolates as sequence with Richardsons method Notes ----- Suppose you have series expansion that goes like this L = f(h) + a0 * h^p_0 + a1 * h^p_1+ a2 * h^p_2 + ... where p_i = order + step * i and f(h) -> L as h -> 0, but f(0) != L. If we evaluate the right hand side for different stepsizes h we can fit a polynomial to that sequence of approximations. This is exactly what this class does. Example ------- >>> import numpy as np >>> import numdifftools as nd >>> n = 3 >>> Ei = np.zeros((n,1)) >>> h = np.zeros((n,1)) >>> linfun = lambda i : np.linspace(0, np.pi/2., 2**(i+5)+1) >>> for k in np.arange(n): ... x = linfun(k) ... h[k] = x[1] ... Ei[k] = np.trapz(np.sin(x),x) >>> En, err, step = nd.Richardson(step=1, order=1)(Ei, h) >>> truErr = Ei-1. >>> (truErr, err, En) (array([[ -2.00805680e-04], [ -5.01999079e-05], [ -1.25498825e-05]]), array([[ 0.00320501]]), array([[ 1.]])) """
[docs] def __init__(self, step_ratio=2.0, step=1, order=1, num_terms=2): self.num_terms = num_terms self.order = order self.step = step self.step_ratio = step_ratio
def _r_matrix(self, num_terms): step = self.step i, j = np.ogrid[0:num_terms + 1, 0:num_terms] r_mat = np.ones((num_terms + 1, num_terms + 1)) r_mat[:, 1:] = (1.0 / self.step_ratio) ** (i * (step * j + self.order)) return r_mat def rule(self, sequence_length=None): if sequence_length is None: sequence_length = self.num_terms + 1 num_terms = min(self.num_terms, sequence_length - 1) if num_terms > 0: r_mat = self._r_matrix(num_terms) return linalg.pinv(r_mat)[0] return np.ones((1,)) @staticmethod def _estimate_error(new_sequence, old_sequence, steps, rule): m, _n = new_sequence.shape if m < 2: return (np.abs(new_sequence) * EPS + steps) * 10.0 cov1 = np.sum(rule**2) # 1 spare dof fact = np.maximum(12.7062047361747 * np.sqrt(cov1), EPS * 10.) err = np.abs(np.diff(new_sequence, axis=0)) * fact tol = np.maximum(np.abs(new_sequence[1:]), np.abs(new_sequence[:-1])) * EPS * fact converged = err <= tol abserr = err + np.where(converged, tol * 10, abs(new_sequence[:-1] - old_sequence[1:]) * fact) return abserr def extrapolate(self, sequence, steps): return self.__call__(sequence, steps) def __call__(self, sequence, steps): ne = sequence.shape[0] rule = self.rule(ne) nr = rule.size - 1 m = ne - nr new_sequence = convolve(sequence, rule[::-1], axis=0, origin=nr // 2) abserr = self._estimate_error(new_sequence, sequence, steps, rule) return new_sequence[:m], abserr[:m], steps[:m]
if __name__ == '__main__': dea_demo() # epsalg_demo()