# -*- coding: utf-8 -*- # type: ignore """Chemical Engineering Design Library (ChEDL). Utilities for process modeling. Copyright (C) 2018, 2019, 2020 Caleb Bell Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicensse, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ from __future__ import division from math import (sin, exp, pi, fabs, copysign, log, isinf, acos, cos, sin, atan2, asinh, sqrt, gamma) from cmath import sqrt as csqrt, log as clog import sys from fluids.numerics.arrays import (solve as py_solve, inv, dot, norm2, inner_product, eye, array_as_tridiagonals, tridiagonals_as_array, solve_tridiagonal, subset_matrix) from fluids.numerics.special import (py_hypot, py_cacos, py_catan, py_catanh, trunc_exp, trunc_log) __all__ = ['isclose', 'horner', 'horner_and_der', 'horner_and_der2', 'horner_and_der3', 'quadratic_from_f_ders', 'chebval', 'interp', 'linspace', 'logspace', 'cumsum', 'diff', 'basic_damping', 'is_poly_negative', 'is_poly_positive', 'implementation_optimize_tck', 'tck_interp2d_linear', 'bisect', 'ridder', 'brenth', 'newton', 'secant', 'halley', 'one_sided_secant', 'splev', 'bisplev', 'derivative', 'jacobian', 'hessian', 'normalize', 'oscillation_checker', 'IS_PYPY', 'roots_cubic', 'roots_quartic', 'newton_system', 'broyden2', 'basic_damping', 'solve_2_direct', 'solve_3_direct', 'solve_4_direct', 'sincos', 'horner_and_der4', 'lambertw', 'ellipe', 'gamma', 'gammaincc', 'erf', 'i1', 'i0', 'k1', 'k0', 'iv', 'mean', 'polylog2', 'roots_quadratic', 'numpy', 'nquad', 'catanh', 'factorial', 'polyint_over_x', 'horner_log', 'polyint', 'zeros', 'full', 'chebder', 'chebint', 'exp_cheb', 'polyder', 'make_damp_initial', 'quadratic_from_points', 'OscillationError', 'UnconvergedError', 'caching_decorator', 'NoSolutionError', 'SamePointError', 'NotBoundedError', 'damping_maintain_sign', 'oscillation_checking_wrapper', 'trunc_exp', 'trunc_log', 'fit_integral_linear_extrapolation', 'fit_integral_over_T_linear_extrapolation', 'poly_fit_integral_value', 'poly_fit_integral_over_T_value', 'evaluate_linear_fits', 'evaluate_linear_fits_d', 'evaluate_linear_fits_d2', 'best_bounding_bounds', 'newton_minimize', 'array_as_tridiagonals', 'tridiagonals_as_array', 'solve_tridiagonal', 'subset_matrix', 'assert_close', 'assert_close1d', 'assert_close2d', 'assert_close3d', 'assert_close4d', 'translate_bound_func', 'translate_bound_jac', 'translate_bound_f_jac', 'curve_fit', 'quad', 'quad_adaptive', 'stable_poly_to_unstable', 'std', 'min_max_ratios', 'detect_outlier_normal', 'max_abs_error', 'max_abs_rel_error', 'max_squared_error', 'max_squared_rel_error', 'mean_abs_error', 'mean_abs_rel_error', 'mean_squared_error', 'mean_squared_rel_error', # Complex number math missing in micropython 'cacos', 'catan', 'deflate_cubic_real_roots', 'fit_minimization_targets', 'root', 'minimize', 'fsolve', 'differential_evolution', 'lmder', 'lmfit', 'horner_backwards', 'exp_horner_backwards', 'horner_backwards_ln_tau', 'exp_horner_backwards_ln_tau', 'exp_horner_backwards_ln_tau_and_der', 'exp_horner_backwards_ln_tau_and_der2', 'exp_poly_ln_tau_coeffs2', 'exp_poly_ln_tau_coeffs3', 'exp_horner_backwards_and_der', 'exp_horner_backwards_and_der2', 'exp_horner_backwards_and_der3', 'horner_backwards_ln_tau_and_der', 'horner_backwards_ln_tau_and_der2', 'horner_backwards_ln_tau_and_der3', 'horner_domain', 'polynomial_offset_scale', 'horner_stable', 'horner_stable_and_der', 'horner_stable_and_der2', 'horner_stable_and_der3', 'horner_stable_and_der4', 'exp_horner_stable', 'exp_horner_stable_and_der', 'exp_horner_stable_and_der2', 'exp_horner_stable_and_der3', 'exp_cheb_and_der', 'exp_cheb_and_der2', 'exp_cheb_and_der3', 'chebval_ln_tau', 'chebval_ln_tau_and_der', 'chebval_ln_tau_and_der2', 'chebval_ln_tau_and_der3', 'horner_stable_ln_tau', 'horner_stable_ln_tau_and_der', 'horner_stable_ln_tau_and_der2', 'horner_stable_ln_tau_and_der3', 'exp_cheb_ln_tau', 'exp_cheb_ln_tau_and_der', 'exp_cheb_ln_tau_and_der2', 'exp_horner_stable_ln_tau', 'exp_horner_stable_ln_tau_and_der', 'exp_horner_stable_ln_tau_and_der2', ] from fluids.numerics import doubledouble from fluids.numerics.doubledouble import * __all__.extend(doubledouble.__all__) __numba_additional_funcs__ = ['py_bisplev', 'py_splev', 'binary_search', 'py_lambertw', '_lambertw_err', 'newton_err', 'norm2', 'py_solve', 'func_35_splev', 'func_40_splev', 'quad_adaptive', 'fixed_quad_Gauss_Kronrod', 'halley_compat_numba', ] nan = float("nan") inf = float("inf") SKIP_DEPENDENCIES = False # for testing class FakePackage(object): pkg = None def __getattr__(self, name): raise ImportError('%s in not installed and required by this feature' %(self.pkg)) def __init__(self, pkg): self.pkg = pkg version_components = sys.version.split('.') PY_MAJOR, PY_MINOR = int(version_components[0]), int(version_components[1]) PY37 = (PY_MAJOR, PY_MINOR) >= (3, 7) try: # The right way imports the platform module which costs to ms to load! # implementation = platform.python_implementation() IS_PYPY = 'PyPy' in sys.version except AttributeError: IS_PYPY = False # Hacks for numba - allow the module to run in different ways try: IS_PYPY = FORCE_PYPY except: pass try: array_if_needed except: array_if_needed = lambda x: x try: is_micropython = sys.implementation.name == 'micropython' if is_micropython: IS_PYPY = True except: is_micropython = False try: is_ironpython = sys.implementation.name == 'ironpython' if is_ironpython: IS_PYPY = True except: is_ironpython = False if is_micropython: hypot = py_hypot def sincos(x): return sin(x), cos(x) try: if IS_PYPY: def sincos(x): # fast implementation based of cephes and go PI4A = 7.85398125648498535156e-1 PI4B = 3.77489470793079817668e-8 PI4C = 2.69515142907905952645e-15 M4PI = 1.273239544735162542821171882678754627704620361328125 #// 4/pi sinSign, cosSign = False, False if x < 0: x = -x sinSign = True j = int(x * M4PI) y = float(j) if j&1 == 1: j += 1 y += 1 j &= 7 if j > 3: j -= 4 sinSign, cosSign = not sinSign, not cosSign if j > 1: cosSign = not cosSign z = ((x - y*PI4A) - y*PI4B) - y*PI4C zz = z * z cos = 1.0 - 0.5*zz + zz*zz*((((((-1.13585365213876817300E-11*zz)+2.08757008419747316778E-9) *zz+-2.75573141792967388112E-7)*zz+2.48015872888517045348E-5) *zz+-1.38888888888730564116E-3)*zz+4.16666666666665929218E-2) sin = z + z*zz*((((((1.58962301576546568060E-10*zz)+-2.50507477628578072866E-8)*zz+2.75573136213857245213E-6) *zz+-1.98412698295895385996E-4)*zz+8.33333333332211858878E-3)*zz+-1.66666666666666307295E-1) if j == 1 or j == 2: sin, cos = cos, sin if cosSign : cos = -cos if sinSign: sin = -sin return sin, cos except: pass try: from cmath import acos as cacos, atan as catan, atanh as catanh, isclose as cisclose except: cacos = py_cacos catan = py_catan catanh = py_catanh _wraps = None def my_wraps(): global _wraps if _wraps is not None: return _wraps from functools import wraps _wraps = wraps return _wraps #IS_PYPY = True # for testing if not SKIP_DEPENDENCIES: try: # Regardless of actual interpreter, fall back to pure python implementations # if scipy and numpy are not available. import numpy np = numpy except ImportError: # Allow a fake numpy to be imported, but will raise an excption on any use numpy = FakePackage('numpy') IS_PYPY = True else: numpy = FakePackage('numpy') IS_PYPY = True IS_PYPY_OR_SKIP_DEPENDENCIES = IS_PYPY or SKIP_DEPENDENCIES np = numpy try: scalar_types = (float, int, np.float64, np.float32, np.float16, np.float128, np.int8, np.int16, np.int32, np.int64) except: scalar_types = (float, int) #IS_PYPY = True try: from sys import float_info epsilon = float_info.epsilon except: # Probably micropython epsilon = 2.220446049250313e-16 one_epsilon_larger = 1.0 + epsilon one_epsilon_smaller = 1.0 - epsilon zero_epsilon_smaller = 0.0 - epsilon one_10_epsilon_larger = 1.0 + epsilon*10.0 one_10_epsilon_smaller = 1.0 - epsilon*10.0 one_100_epsilon_larger = 1.0 + epsilon*100.0 one_100_epsilon_smaller = 1.0 - epsilon*100.0 one_1000_epsilon_larger = 1.0 + epsilon*1000.0 one_1000_epsilon_smaller = 1.0 - epsilon*1000.0 _iter = 100 _xtol = 1e-12 _rtol = epsilon*2.0 third = 1.0/3.0 sixth = 1.0/6.0 ninth = 1.0/9.0 twelfth = 1.0/12.0 two_thirds = 2.0/3.0 four_thirds = 4.0/3.0 root_three = 1.7320508075688772 # sqrt(3.0) one_27 = 1.0/27.0 complex_factor = 0.8660254037844386j # (sqrt(3)*0.5j) def roots_quadratic(a, b, c): if a == 0.0: return (-c/b, ) D = b*b - 4.0*a*c a_inv_2 = 0.5/a if D < 0.0: D = sqrt(-D) x1 = (-b + D*1.0j)*a_inv_2 x2 = (-b - D*1.0j)*a_inv_2 else: D = sqrt(D) x1 = (D - b)*a_inv_2 x2 = -(b + D)*a_inv_2 return (x1, x2) def roots_cubic_a1(b, c, d): # Output from mathematica t1 = b*b t2 = t1*b t4 = c*b t9 = c*c t16 = d*d t19 = csqrt(12.0*t9*c + 12.0*t2*d - 54.0*t4*d - 3.0*t1*t9 + 81.0*t16) t22 = (-8.0*t2 + 36.0*t4 - 108.0*d + 12.0*t19)**third root1 = t22*sixth - 6.0*(c*third - t1*ninth)/t22 - b*third t28 = (c*third - t1*ninth)/t22 t101 = -t22*twelfth + 3.0*t28 - b*third t102 = root_three*(t22*sixth + 6.0*t28) root2 = t101 + 0.5j*t102 root3 = t101 - 0.5j*t102 return [root1, root2, root3] def roots_cubic_a2(a, b, c, d): # Output from maple t2 = a*a t3 = d*d t10 = c*c t14 = b*b t15 = t14*b t20 = csqrt(-18.0*a*b*c*d + 4.0*a*t10*c + 4.0*t15*d - t14*t10 + 27.0*t2*t3) t31 = (36.0*c*b*a + 12.0*root_three*t20*a - 108.0*d*t2 - 8.0*t15)**third t32 = 1.0/a root1 = t31*t32*sixth - two_thirds*(3.0*a*c - t14)*t32/t31 - b*t32*third t33 = t31*t32 t40 = (3.0*a*c - t14)*t32/t31 t50 = -t33*twelfth + t40*third - b*t32*third t51 = 0.5j*root_three *(t33*sixth + two_thirds*t40) root2 = t50 + t51 root3 = t50 - t51 return [root1, root2, root3] def roots_cubic(a, b, c, d): r'''Cubic equation solver based on a variety of sources, algorithms, and numerical tools. It seems evident after some work that no analytical solution using floating points will ever have high-precision results for all cases of inputs. Some other solvers, such as NumPy's roots which uses eigenvalues derived using some BLAS, seem to provide bang-on answers for all values coefficients. However, they can be quite slow - and where possible there is still a need for analytical solutions to obtain 15-35x speed, such as when using PyPy. A particular focus of this routine is where a=1, b is a small number in the range -10 to 10 - a common occurrence in cubic equations of state. Parameters ---------- a : float Coefficient of x^3, [-] b : float Coefficient of x^2, [-] c : float Coefficient of x, [-] d : float Added coefficient, [-] Returns ------- roots : tuple(float) The evaluated roots of the polynomial, 1 value when a and b are zero, two values when a is zero, and three otherwise, [-] Notes ----- For maximum speed, provide Python floats. Compare the speed with numpy via: %timeit roots_cubic(1.0, 100.0, 1000.0, 10.0) %timeit np.roots([1.0, 100.0, 1000.0, 10.0]) %timeit roots_cubic(1.0, 2.0, 3.0, 4.0) %timeit np.roots([1.0, 2.0, 3.0, 4.0]) The speed is ~15-35 times faster; or using PyPy, 240-370 times faster. Examples -------- >>> roots_cubic(1.0, 100.0, 1000.0, 10.0) (-0.0100100190, -88.731288, -11.25870159) References ---------- .. [1] "Solving Cubic Equations." Accessed January 5, 2019. http://www.1728.org/cubic2.htm. ''' ''' Notes ----- Known issue is inputs that look like 1, -0.999999999978168, 1.698247818501352e-11, -8.47396642608142e-17 Errors grown unbound, starting when b is -.99999 and close to 1. ''' if a == 0.0: if b == 0.0: return (-d/c, ) D = c*c - 4.0*b*d b_inv_2 = 0.5/b if D < 0.0: D = sqrt(-D) x1 = (-c + D*1.0j)*b_inv_2 x2 = (-c - D*1.0j)*b_inv_2 else: D = sqrt(D) x1 = (D - c)*b_inv_2 x2 = -(c + D)*b_inv_2 return (x1, x2) a_inv = 1.0/a a_inv2 = a_inv*a_inv bb = b*b '''Herbie modifications for f: c*a_inv - b_a*b_a*third ''' b_a = b*a_inv b_a2 = b_a*b_a f = c*a_inv - b_a2*third # f = (3.0*c*a_inv - bb*a_inv2)*third g = ((2.0*(bb*b) * a_inv2*a_inv) - (9.0*b*c)*(a_inv2) + (27.0*d*a_inv))*one_27 # g = (((2.0/(a/b))/((a/b) * (a/b)) + d*27.0/a) - (9.0/a*b)*c/a)/27.0 h = (0.25*(g*g) + (f*f*f)*one_27) # print(f, g, h) '''h has no savings on precision - 0.4 error to 0.2. ''' # print(f, g, h, 'f, g, h') if h == 0.0 and g == 0.0 and f == 0.0: if d/a >= 0.0: x = -((d*a_inv)**(third)) else: x = (-d*a_inv)**(third) return (x, x, x) elif h > 0.0: # Happy with these formulas - double doubles should be fast. # No complex numbers are needed here. # print('basic') # 1 real root, 2 imag root_h = sqrt(h) R = -0.5*g + root_h # It is possible to save one of the power of thirds! if R >= 0.0: S = R**third else: S = -((-R)**third) T = -(0.5*g) - root_h if T >= 0.0: U = (T**(third)) else: U = -(((-T)**(third))) SU = S + U b_3a = b*(third*a_inv) t1 = -0.5*SU - b_3a t2 = (S - U)*complex_factor x1 = SU - b_3a # x1 is OK actually in some tests? the issue is x2, x3? x2 = t1 + t2 x3 = t1 - t2 else: # elif h <= 0.0: t2 = a*a t3 = d*d t10 = c*c t14 = b*b t15 = t14*b '''This method is inaccurate when choice_term is too small; but still more accurate than the other method. ''' choice_term = -18.0*a*b*c*d + 4.0*a*t10*c + 4.0*t15*d - t14*t10 + 27.0*t2*t3 if (abs(choice_term) > 1e-12 or abs(b + 1.0) < 1e-7): # print('mine') t32 = 1.0/a t20 = csqrt(choice_term) t31 = (36.0*c*b*a + 12.0*root_three*t20*a - 108.0*d*t2 - 8.0*t15)**third t33 = t31*t32 t32_t31 = t32/t31 x1 = (t33*sixth - two_thirds*(3.0*a*c - t14)*t32_t31 - b*t32*third).real t40 = (3.0*a*c - t14)*t32_t31 t50 = -t33*twelfth + t40*third - b*t32*third t51 = 0.5j*root_three*(t33*sixth + two_thirds*t40) x2 = (t50 + t51).real x3 = (t50 - t51).real else: # print('other') # 3 real roots # example is going in here i = sqrt(((g*g)*0.25) - h) j = i**third # There was a saving for j but it was very weird with if statements! '''Clamied nothing saved for k. ''' k = acos(-0.5*g/i) # L = -j # N, M = sincos(k*third) # N *= root_three k_third = k*third M = cos(k_third) N = root_three*sin(k_third) P = -b_a*third # Direct formula for x1 x1 = 2.0*j*M + P x2 = P - j*(M + N) x3 = P - j*(M - N) return (x1, x2, x3) def roots_quartic(a, b, c, d, e): # There is no divide by zero check. A should be 1 for best numerical results # Multiple order of magnitude differences still can cause problems # Like [1, 0.0016525874561771799, 106.8665062954208, 0.0032802613917246727, 0.16036091315844248] x0 = 1.0/a x1 = b*x0 x2 = -x1*0.25 x3 = c*x0 x4 = b*b*x0*x0 x5 = -two_thirds*x3 + 0.25*x4 x6 = x3 - 0.375*x4 x6_2 = x6*x6 x7 = x6_2*x6 x8 = d*x0 x9 = x1*(-0.5*x3 + 0.125*x4) x10 = (x8 + x9)*(x8 + x9) x11 = e*x0 x12 = x1*(x1*(-0.0625*x3 + 0.01171875*x4) + 0.25*x8) # 0.01171875 = 3/256 x13 = x6*(x11 - x12) x14 = -.125*x10 + x13*third - x7/108.0 x15 = 2.0*(x14 + 0.0j)**(third) x16 = csqrt(-x15 + x5) x17 = 0.5*x16 x18 = -x17 + x2 x19 = -four_thirds*x3 x20 = 0.5*x4 x21 = x15 + x19 + x20 x22 = 2.0*x8 + 2.0*x9 x23 = x22/x16 x24 = csqrt(x21 + x23)*0.5 x25 = -x11 + x12 - twelfth*x6_2 x27 = (0.0625*x10 - x13*sixth + x7/216.0 + csqrt(0.25*x14*x14 + one_27*x25*x25*x25))**(third) x28 = 2.0*x27 x29 = two_thirds*x25/(x27) x30 = csqrt(x28 - x29 + x5) x31 = 0.5*x30 x32 = x2 - x31 x33 = x19 + x20 - x28 + x29 x34 = x22/x30 x35 = csqrt(x33 + x34)*0.5 x36 = x17 + x2 x37 = csqrt(x21 - x23)*0.5 x38 = x2 + x31 x39 = csqrt(x33 - x34)*0.5 return ((x32 - x35), (x32 + x35), (x38 - x39), (x38 + x39)) def mean(data): # Much faster than the statistics.mean module return sum(data)/len(data) def detect_outlier_normal(data, cutoff=3): std_val = std(data) mean_val = mean(data) low = mean_val - std_val*cutoff high = mean_val + std_val*cutoff outlier_idxs = [] for i, p in enumerate(data): if p < low or p > high: outlier_idxs.append(i) return outlier_idxs def linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None): """Port of numpy's linspace to pure python. Does not support dtype, and returns lists of floats. """ num = int(num) start = start * 1. stop = stop * 1. if num <= 0: return [] if endpoint: if num == 1: return [start] step = (stop-start)/float((num-1)) if num == 1: step = nan y = [start] for _ in range(num-2): y.append(y[-1] + step) y.append(stop) else: step = (stop-start)/float(num) if num == 1: step = nan y = [start] for _ in range(num-1): y.append(y[-1] + step) if retstep: return y, step else: return y def logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None): y = linspace(start, stop, num=num, endpoint=endpoint) for i in range(len(y)): y[i] = base**y[i] # return [base**yi for yi in y] return y def product(l): # Helper in some functions tot = 1.0 for i in l: tot *= i return tot def cumsum(a): # Does not support multiple dimensions base = a[0] sums = [base] for v in a[1:]: base = base + v sums.append(base) return sums def diff(a, n=1, axis=-1): if n == 0: return a if n < 0: raise ValueError( "order must be non-negative but got %s" %(n)) # nd = 1 # hardcode diffs = [] for i in range(1, len(a)): delta = a[i] - a[i-1] diffs.append(delta) if n > 1: return diff(diffs, n-1) return diffs central_diff_weights_precomputed = { (1, 3): [-0.5, 0.0, 0.5], (1, 5): [0.08333333333333333, -0.6666666666666666, 0.0, 0.6666666666666666, -0.08333333333333333], (1, 7): [-0.016666666666666666, 0.15, -0.75, 0.0, 0.75, -0.15, 0.016666666666666666], (1, 9): [0.0035714285714285713, -0.0380952380952381, 0.2, -0.8, 0.0, 0.8, -0.2, 0.0380952380952381, -0.0035714285714285713], (1, 11): [-0.0007936507936507937, 0.00992063492063492, -0.05952380952380952, 0.23809523809523808, -0.8333333333333334, 0.0, 0.8333333333333334, -0.23809523809523808, 0.05952380952380952, -0.00992063492063492, 0.0007936507936507937], (1, 13): [0.00018037518037518038, -0.0025974025974025974, 0.017857142857142856, -0.07936507936507936, 0.26785714285714285, -0.8571428571428571, 0.0, 0.8571428571428571, -0.26785714285714285, 0.07936507936507936, -0.017857142857142856, 0.0025974025974025974, -0.00018037518037518038], (1, 15): [-4.1625041625041625e-05, 0.0006798756798756799, -0.005303030303030303, 0.026515151515151516, -0.09722222222222222, 0.2916666666666667, -0.875, 0.0, 0.875, -0.2916666666666667, 0.09722222222222222, -0.026515151515151516, 0.005303030303030303, -0.0006798756798756799, 4.1625041625041625e-05], (1, 17): [9.712509712509713e-06, -0.0001776001776001776, 0.001554001554001554, -0.008702408702408702, 0.03535353535353535, -0.11313131313131314, 0.3111111111111111, -0.8888888888888888, 0.0, 0.8888888888888888, -0.3111111111111111, 0.11313131313131314, -0.03535353535353535, 0.008702408702408702, -0.001554001554001554, 0.0001776001776001776, -9.712509712509713e-06], (1, 19): [-2.285296402943462e-06, 4.6277252159605104e-05, -0.00044955044955044955, 0.002797202797202797, -0.012587412587412588, 0.044055944055944055, -0.12727272727272726, 0.32727272727272727, -0.9, 0.0, 0.9, -0.32727272727272727, 0.12727272727272726, -0.044055944055944055, 0.012587412587412588, -0.002797202797202797, 0.00044955044955044955, -4.6277252159605104e-05, 2.285296402943462e-06], (2, 3): [1.0, -2.0, 1.0], (2, 5): [-0.08333333333333333, 1.3333333333333333, -2.5, 1.3333333333333333, -0.08333333333333333], (2, 7): [0.011111111111111112, -0.15, 1.5, -2.7222222222222223, 1.5, -0.15, 0.011111111111111112], (2, 9): [-0.0017857142857142857, 0.025396825396825397, -0.2, 1.6, -2.8472222222222223, 1.6, -0.2, 0.025396825396825397, -0.0017857142857142857], (2, 11): [0.00031746031746031746, -0.00496031746031746, 0.03968253968253968, -0.23809523809523808, 1.6666666666666667, -2.9272222222222224, 1.6666666666666667, -0.23809523809523808, 0.03968253968253968, -0.00496031746031746, 0.00031746031746031746], (2, 13): [-6.012506012506013e-05, 0.001038961038961039, -0.008928571428571428, 0.05291005291005291, -0.26785714285714285, 1.7142857142857142, -2.9827777777777778, 1.7142857142857142, -0.26785714285714285, 0.05291005291005291, -0.008928571428571428, 0.001038961038961039, -6.012506012506013e-05], (2, 15): [1.1892869035726179e-05, -0.00022662522662522663, 0.0021212121212121214, -0.013257575757575758, 0.06481481481481481, -0.2916666666666667, 1.75, -3.02359410430839, 1.75, -0.2916666666666667, 0.06481481481481481, -0.013257575757575758, 0.0021212121212121214, -0.00022662522662522663, 1.1892869035726179e-05], (2, 17): [-2.428127428127428e-06, 5.074290788576503e-05, -0.000518000518000518, 0.003480963480963481, -0.017676767676767676, 0.07542087542087542, -0.3111111111111111, 1.7777777777777777, -3.05484410430839, 1.7777777777777777, -0.3111111111111111, 0.07542087542087542, -0.017676767676767676, 0.003480963480963481, -0.000518000518000518, 5.074290788576503e-05, -2.428127428127428e-06], (2, 19): [5.078436450985471e-07, -1.1569313039901276e-05, 0.00012844298558584272, -0.0009324009324009324, 0.005034965034965035, -0.022027972027972027, 0.08484848484848485, -0.32727272727272727, 1.8, -3.0795354623330815, 1.8, -0.32727272727272727, 0.08484848484848485, -0.022027972027972027, 0.005034965034965035, -0.0009324009324009324, 0.00012844298558584272, -1.1569313039901276e-05, 5.078436450985471e-07], (3, 5): [-0.5, 1.0, 0.0, -1.0, 0.5], (3, 7): [0.125, -1.0, 1.625, 0.0, -1.625, 1.0, -0.125], (3, 9): [-0.029166666666666667, 0.30000000000000004, -1.4083333333333332, 2.033333333333333, 0.0, -2.033333333333333, 1.4083333333333332, -0.30000000000000004, 0.029166666666666667], # (3, 11): [0.006779100529100529, -0.08339947089947089, 0.48303571428571435, -1.7337301587301588, 2.3180555555555555, # 0.0, -2.3180555555555555, 1.7337301587301588, -0.48303571428571435, 0.08339947089947089, # -0.006779100529100529], # (3, 13): [-0.0015839947089947089, 0.02261904761904762, -0.1530952380952381, 0.6572751322751322, -1.9950892857142857, # 2.527142857142857, 0.0, -2.527142857142857, 1.9950892857142857, -0.6572751322751322, 0.1530952380952381, # -0.02261904761904762, 0.0015839947089947089], # (3, 15): [0.0003724747474747475, -0.006053691678691679, 0.04682990620490621, -0.2305699855699856, 0.8170667989417989, # -2.2081448412698412, 2.6869345238095237, 0.0, -2.6869345238095237, 2.2081448412698412, -0.8170667989417989, # 0.2305699855699856, -0.04682990620490621, 0.006053691678691679, -0.0003724747474747475], # (3, 17): [-8.810006131434702e-05, 0.0016058756058756059, -0.013982697196982911, 0.07766492766492766, # -0.31074104136604136, 0.9613746993746994, -2.384521164021164, 2.81291761148904, 0.0, -2.81291761148904, # 2.384521164021164, -0.9613746993746994, 0.31074104136604136, -0.07766492766492766, 0.013982697196982911, # -0.0016058756058756059, 8.810006131434702e-05], # (3, 19): [2.0943672729387014e-05, -0.00042319882498453924, 0.00409817266067266, -0.025376055161769447, # 0.11326917130488559, -0.3904945471195471, 1.0909741462241462, -2.532634817563389, 2.9147457482993198, 0.0, # -2.9147457482993198, 2.532634817563389, -1.0909741462241462, 0.3904945471195471, -0.11326917130488559, # 0.025376055161769447, -0.00409817266067266, 0.00042319882498453924, -2.0943672729387014e-05], # (4, 5): [1.0, -4.0, 6.0, -4.0, 1.0], # (4, 7): [-0.16666666666666666, 2.0, -6.5, 9.333333333333334, -6.5, 2.0, -0.16666666666666666], # (4, 9): [0.029166666666666667, -0.4, 2.8166666666666664, -8.133333333333333, 11.375, -8.133333333333333, # 2.8166666666666664, -0.4, 0.029166666666666667], # (4, 11): [-0.005423280423280424, 0.08339947089947089, -0.644047619047619, 3.4674603174603176, -9.272222222222222, # 12.741666666666665, -9.272222222222222, 3.4674603174603176, -0.644047619047619, 0.08339947089947089, # -0.005423280423280424], # (4, 13): [0.0010559964726631393, -0.018095238095238095, 0.1530952380952381, -0.8763668430335096, 3.9901785714285714, # -10.108571428571429, 13.717407407407407, -10.108571428571429, 3.9901785714285714, -0.8763668430335096, # 0.1530952380952381, -0.018095238095238095, 0.0010559964726631393], # (5, 7): [-0.5, 2.0, -2.5, 0.0, 2.5, -2.0, 0.5], # (5, 9): [0.16666666666666669, -1.5, 4.333333333333333, -4.833333333333334, 0.0, 4.833333333333334, -4.333333333333333, # 1.5, -0.16666666666666669], # (5, 11): [-0.04513888888888889, 0.5277777777777778, -2.71875, 6.5, -6.729166666666667, 0.0, 6.729166666666667, -6.5, # 2.71875, -0.5277777777777778, 0.04513888888888889], # (5, 13): [0.011491402116402117, -0.16005291005291006, 1.033399470899471, -3.9828042328042326, 8.39608134920635, # -8.246031746031747, 0.0, 8.246031746031747, -8.39608134920635, 3.9828042328042326, -1.033399470899471, # 0.16005291005291006, -0.011491402116402117] } def deflate_cubic_real_roots(b, c, d, x0): F = b + x0 G = -d/x0 D = F*F - 4.0*G # if D < 0.0: # D = (-D)**0.5 # x1 = (-F + D*1.0j)*0.5 # x2 = (-F - D*1.0j)*0.5 # else: if D < 0.0: return (0.0, 0.0) D = sqrt(D) x1 = 0.5*(D - F)#(D - c)*0.5 x2 = 0.5*(-F - D) #-(c + D)*0.5 return x1, x2 def central_diff_weights(points, divisions=1): # Check the cache if (divisions, points) in central_diff_weights_precomputed: return central_diff_weights_precomputed[(divisions, points)] if points < divisions + 1: raise ValueError("Points < divisions + 1, cannot compute") if points % 2 == 0: raise ValueError("Odd number of points required") ho = points >> 1 x = [[xi] for xi in range(-ho, ho+1)] X = [] for xi in x: line = [1.0] + [xi[0]**k for k in range(1, points)] X.append(line) factor = product(range(1, divisions + 1)) # from scipy.linalg import inv from sympy import Matrix # The above coefficients were generated from a routine which used Fractions # Additional coefficients cannot reliable be computed with numpy or floating # point numbers - the error is too great # sympy must be used for reliability inverted = [[float(j) for j in i] for i in Matrix(X).inv().tolist()] w = [i*factor for i in inverted[divisions]] # w = [i*factor for i in (inv(X)[divisions]).tolist()] central_diff_weights_precomputed[(divisions, points)] = w return w def derivative(func, x0, dx=1.0, n=1, args=(), order=3, scalar=True, lower_limit=None, upper_limit=None, kwargs={}): """Reimplementation of SciPy's derivative function, with more cached coefficients and without using numpy. If new coefficients not cached are needed, they are only calculated once and are remembered. Support for vector value functions has also been added. """ if order < n + 1: raise ValueError if order % 2 == 0: raise ValueError weights = central_diff_weights(order, n) ho = order >> 1 denominator = 1.0/product([dx]*n) if scalar: max_x = x0 + (order - 1 - ho)*dx if upper_limit is not None and max_x > upper_limit: x0 -= (max_x - x0) min_x = x0 + -ho*dx if lower_limit is not None and min_x < lower_limit: x0 += (x0 - min_x) tot = 0.0 for k in range(order): if weights[k] != 0.0: tot += weights[k]*func(x0 + (k - ho)*dx, *args, **kwargs) return tot*denominator else: numerators = None for k in range(order): f = func(x0 + (k - ho)*dx, *args, **kwargs) if numerators is None: N = len(f) numerators = [0.0]*N for i in range(N): numerators[i] += weights[k]*f[i] return [num*denominator for num in numerators] kronrod_weights = { # 2: [0.1979797979797974, 0.4909090909090922, 0.6222222222222221, 0.4909090909090904, 0.1979797979797974], # 3: [0.10465622602646653, 0.2684880898683364, 0.4013974147759613, 0.45091653865847436, 0.40139741477596186, 0.26848808986833306, 0.10465622602646701], # 4: [0.06297737366547282, 0.17005360533572325, 0.26679834045228396, 0.3269491896014517, 0.3464429818901376, 0.32694918960145064, 0.2667983404522837, # 0.17005360533572333, 0.06297737366547272 # ], # 5: [0.042582036751078814, 0.11523331662247298, 0.1868007965564952, 0.2410403392286488, 0.27284980191255903, 0.28298741785749215, 0.2728498019125575, # 0.24104033922864848, 0.1868007965564929, 0.11523331662247287, 0.04258203675108158 # ], # 6: [0.03039615411981852, 0.08369444044690531, 0.1373206046344473, 0.1810719943231391, 0.21320965227196226, 0.23377086411699666, 0.2410725801734652, # 0.23377086411699363, 0.21320965227196081, 0.1810719943231382, 0.13732060463444687, 0.08369444044690623, 0.030396154119819826 # ], # 7: [0.02293532201053038, 0.0630920926299766, 0.10479001032225017, 0.14065325971552478, 0.16900472663927035, 0.19035057806478545, 0.2044329400753001, # 0.20948214108472712, 0.2044329400752992, 0.1903505780647853, 0.16900472663926686, 0.14065325971552592, 0.10479001032225005, 0.06309209262997834, # 0.022935322010529363 # ], # 8: [0.017822383320711625, 0.04943939500213785, 0.08248229893135677, 0.11164637082684088, 0.13626310925517557, 0.15665260616818655, 0.17207060855521186, # 0.18140002506803385, 0.18444640574468968, 0.1814000250680349, 0.1720706085552107, 0.1566526061681889, 0.13626310925517285, 0.11164637082683933, # 0.08248229893135825, 0.04943939500213947, 0.017822383320710476 # ], # 9: [0.014304775643840032, 0.03963189516026222, 0.0665181559402729, 0.09079068168872734, 0.11178913468441762, 0.13000140685534053, 0.14523958838436735, # 0.15641352778848416, 0.16286282744011463, 0.16489601282834765, 0.16286282744011538, 0.15641352778848497, 0.14523958838436518, # 0.13000140685534053, 0.11178913468441844, 0.09079068168872709, 0.0665181559402743, 0.03963189516026152, 0.01430477564383906 # ], 10: [0.011694638867371172, 0.03255816230796458, 0.05475589657435185, 0.07503967481091954, 0.09312545458369915, 0.10938715880229993, 0.12349197626206619, 0.1347092173114735, 0.1427759385770589, 0.14773910490133732, 0.1494455540029151, 0.14773910490133846, 0.14277593857705978, 0.13470921731147467, 0.1234919762620664, 0.10938715880229763, 0.09312545458369788, 0.0750396748109198, 0.05475589657435227, 0.03255816230796462, 0.01169463886737192 ], # 11: [0.009765441045959541, 0.02715655468210335, 0.045829378564428196, 0.06309742475037379, 0.07866457193222803, 0.092953098596902, 0.1058720744813897, # 0.11673950246104702, 0.12515879910031777, 0.13128068422980624, 0.13519357279988428, 0.1365777947111174, 0.13519357279988467, 0.1312806842298052, # 0.12515879910031882, 0.11673950246104742, 0.10587207448139022, 0.09295309859690085, 0.07866457193222773, 0.06309742475037532, # 0.04582937856442653, 0.027156554682104143, 0.009765441045960747 # ], # 12: [0.008257711433166826, 0.023036084038981483, 0.03891523046929862, 0.053697017607756393, 0.06725090705084145, 0.07992027533360162, # 0.09154946829504902, 0.1016497322790617, 0.11002260497764366, 0.11671205350175794, 0.12162630352394828, 0.12458416453615645, # 0.12555689390547403, 0.12458416453615662, 0.12162630352394885, 0.11671205350175637, 0.11002260497764366, 0.10164973227906021, # 0.09154946829504902, 0.07992027533360205, 0.06725090705084015, 0.05369701760775601, 0.038915230469299136, 0.023036084038982187, # 0.008257711433168356 # ], # 13: [0.007087846351247695, 0.019753746382707126, 0.03344358998955114, 0.046279017973829016, 0.05811521042311424, 0.06930363324778273, # 0.0798059621694777, 0.08916844187754189, 0.09714173487607722, 0.10383060116904075, 0.10926635109528536, 0.11321025917152877, 0.1154887990912863, # 0.11620961236306047, 0.11548879909128668, 0.11321025917152944, 0.10926635109528575, 0.10383060116903994, 0.09714173487607768, # 0.08916844187754075, 0.0798059621694761, 0.06930363324778134, 0.05811521042311445, 0.0462790179738303, 0.033443589989552346, # 0.019753746382705984, 0.007087846351248617 # ], # 14: [0.006139558686377161, 0.01714845890993743, 0.02904870126150814, 0.04025059487268863, 0.050691543260464045, 0.06066712586674275, # 0.07010297900274666, 0.0786557972496211, 0.08618376286946026, 0.09273683001785471, 0.09826492647210357, 0.1026166273213988, 0.10573163984176341, # 0.10762642111411824, 0.10827006650642954, 0.10762642111411819, 0.10573163984176337, 0.10261662732139956, 0.09826492647210401, # 0.09273683001785424, 0.08618376286945752, 0.07865579724962145, 0.07010297900274708, 0.060667125866742444, 0.050691543260465384, # 0.040250594872688804, 0.029048701261508613, 0.01714845890993556, 0.00613955868637814 # ], # 15: [0.005377479872922721, 0.015007947329316597, 0.02546084732671451, 0.03534636079137684, 0.044589751324763977, 0.053481524690928775, # 0.06200956780067099, 0.06985412131872938, 0.07684968075772279, 0.08308050282313152, 0.08856444305621097, 0.09312659817082411, # 0.09664272698362378, 0.09917359872179198, 0.10076984552387584, 0.10133000701479176, 0.10076984552387498, 0.09917359872179111, # 0.09664272698362311, 0.09312659817082504, 0.08856444305621132, 0.08308050282313374, 0.07684968075772115, 0.06985412131872824, # 0.06200956780067089, 0.05348152469092814, 0.04458975132476507, 0.03534636079137597, 0.0254608473267153, 0.01500794732931606, # 0.005377479872923303 # ], # 16: [0.0047427770492463476, 0.01325793068809046, 0.022498859440048875, 0.031260543647380616, 0.03951295120242179, 0.04750621597640688, # 0.055205633095423194, 0.062358806011834494, 0.06886299519153054, 0.07476982388560004, 0.08005394126371923, 0.08459580379259089, # 0.08833750257911316, 0.0912920328281922, 0.09343867406092181, 0.09472840124723077, 0.09515421608049866, 0.09472840124722949, # 0.09343867406092161, 0.09129203282819215, 0.08833750257911213, 0.0845958037925905, 0.08005394126371795, 0.07476982388559944, # 0.06886299519153126, 0.062358806011835174, 0.05520563309542224, 0.047506215976407244, 0.03951295120242178, 0.031260543647380366, # 0.02249885944004935, 0.013257930688091162, 0.0047427770492472505 # ], # 17: [0.00421897579377578, 0.011785837562288765, 0.020022233953294232, 0.027856722457863144, 0.035249746751878815, 0.04244263020500057, # 0.04943614182396047, 0.055994044530093004, 0.06200091526823072, 0.06752801671813287, 0.07258989011419156, 0.07705622304620212, # 0.08083667844081736, 0.08397257199123528, 0.08649053220565191, 0.08830878229760371, 0.08936184586778975, 0.08969642194398147, # 0.08936184586778813, 0.08830878229760494, 0.08649053220565157, 0.0839725719912349, 0.08083667844081745, 0.07705622304620285, # 0.07258989011419015, 0.06752801671813129, 0.06200091526822982, 0.05599404453009317, 0.04943614182395916, 0.042442630205000643, # 0.03524974675187985, 0.027856722457863275, 0.020022233953294787, 0.011785837562289191, 0.004218975793776932 # ], # 18: [0.00377351012843988, 0.010554430567545547, 0.01793337620570976, 0.024964402753740674, 0.03163401984884463, 0.0381537583251753, # 0.044509195165090054, 0.0505075715021042, 0.056072390382993734, 0.061253565791358926, 0.06603995819281894, 0.07033745002422243, # 0.07409702738131525, 0.0773320701646388, 0.08003001949509578, 0.08214399978333053, 0.08365485143716135, 0.08456933354407732, # 0.08487813861267707, 0.08456933354407585, 0.08365485143716228, 0.08214399978332948, 0.08003001949509637, 0.07733207016463836, # 0.07409702738131455, 0.07033745002422402, 0.06603995819281755, 0.06125356579135898, 0.05607239038299394, 0.05050757150210246, # 0.044509195165090026, 0.03815375832517691, 0.031634019848844154, 0.024964402753740036, 0.017933376205711227, 0.01055443056754483, # 0.0037735101284401234 # ], # 19: [0.0033981531196335965, 0.009499205461711157, 0.016153414369089493, 0.02250869575220442, 0.028544238117114772, 0.03446235286288065, # 0.04027002324951411, 0.04578597308792263, 0.05092622550624418, 0.05575220103348164, 0.060277411503209075, 0.0644023522619292, # 0.06805978049813398, 0.07128142476493002, 0.07408407453373855, 0.07640286193981659, 0.07818691068574178, 0.07946568358535473, # 0.08026621602085142, 0.08054560329299826, 0.08026621602085034, 0.07946568358535418, 0.07818691068574216, 0.07640286193981571, # 0.07408407453373904, 0.07128142476492985, 0.06805978049813312, 0.06440235226192832, 0.06027741150320817, 0.055752201033481275, # 0.05092622550624404, 0.04578597308792251, 0.04027002324951481, 0.03446235286288103, 0.028544238117114873, 0.02250869575220419, # 0.016153414369091356, 0.009499205461710144, 0.0033981531196347167 # ], # 87: [0.00016881480210456692, 0.00047314868570756235, 0.0008076859860306628, 0.0011322615841296547, 0.00144806597820222, 0.0017668191481724435, 0.0020903160499537304, 0.002411576726915393, 0.0027285063230227296, 0.0030454007869937876, 0.0033639128646212944, 0.0036808253607571195, 0.003994559817314871, 0.004307345076427869, 0.004620335406184352, 0.004931600306915534, 0.00523996592505234, 0.005546814695915935, 0.00585298027614558, 0.006157129826902876, 0.006458347449208747, 0.0067575841722384, 0.007055469839437327, 0.007351002065619594, 0.007643435225722091, 0.007933469441304792, 0.008221599308725454, 0.008507026852701362, 0.008789124673244655, 0.009068433903060025, 0.00934535574073612, 0.00961922899522448, 0.009889512165959614, 0.010156639610742505, 0.01042094604150296, 0.010681867318519509, 0.010938926800718108, 0.011192484212786687, 0.011442825936038695, 0.011689459907180238, 0.01193196007256122, 0.012170632572570447, 0.01240572822845446, 0.012636810511622156, 0.012863493923529042, 0.013086045502712524, 0.013304689821902624, 0.013519034459386238, 0.013728727196381123, 0.013934006347263819, 0.014135077241799326, 0.0143315833864912, 0.01452320041557479, 0.01471014561481331, 0.014892610489607272, 0.015070268437801156, 0.01524281878982064, 0.015410463680651272, 0.01557338508816692, 0.015731281717590702, 0.01588387332676284, 0.016031351508057137, 0.01617389219804666, 0.016311215832305204, 0.016443059955062258, 0.016569609315104587, 0.016691036698090543, 0.016807081400707274, 0.01691749655723537, 0.017022463108090843, 0.017122153135412682, 0.017216322428248873, 0.017304737826468557, 0.017387579005841453, 0.0174650194618558, 0.017536829459282623, 0.017602787871331006, 0.01766307528646667, 0.017717868485658792, 0.017766950435146724, 0.017810110503721263, 0.01784753209761472, 0.017879396979083733, 0.017905499199739466, 0.017925637157192337, 0.017939998782076135, 0.017948772394930045, 0.017951761607250333, 0.01794877239493056, 0.017939998782076676, 0.01792563715719179, 0.017905499199739282, 0.017879396979084167, 0.017847532097614346, 0.01781011050372197, 0.017766950435146787, 0.017717868485658865, 0.01766307528646623, 0.017602787871330267, 0.017536829459281828, 0.01746501946185603, 0.01738757900584135, 0.017304737826468315, 0.017216322428249966, 0.01712215313541374, 0.017022463108090982, 0.016917496557235175, 0.016807081400708086, 0.016691036698090224, 0.016569609315104577, 0.016443059955061293, 0.016311215832305943, 0.016173892198046793, 0.01603135150805684, 0.015883873326762922, 0.0157312817175899, 0.015573385088166592, 0.015410463680651668, 0.015242818789820901, 0.015070268437801513, 0.01489261048960673, 0.014710145614813493, 0.014523200415575395, 0.014331583386491533, 0.01413507724179865, 0.013934006347263444, 0.013728727196380418, 0.013519034459385928, 0.01330468982190219, 0.013086045502712456, 0.01286349392352816, 0.012636810511621802, 0.012405728228453924, 0.012170632572569789, 0.011931960072561447, 0.011689459907180488, 0.011442825936038951, 0.011192484212786695, 0.010938926800718342, 0.010681867318518458, 0.010420946041502708, 0.01015663961074195, 0.009889512165958998, 0.009619228995224303, 0.009345355740736385, 0.009068433903060022, 0.00878912467324535, 0.008507026852701478, 0.008221599308725935, 0.007933469441304955, 0.00764343522572274, 0.007351002065619289, 0.007055469839437192, 0.006757584172238039, 0.006458347449209268, 0.0061571298269029575, 0.005852980276145424, 0.00554681469591638, 0.005239965925052114, 0.004931600306916382, 0.00462033540618409, 0.004307345076429006, 0.003994559817315303, 0.0036808253607556987, 0.00336391286462138, 0.0030454007869938063, 0.0027285063230224034, 0.0024115767269154683, 0.0020903160499540895, 0.0017668191481723982, 0.0014480659782023144, 0.0011322615841305992, 0.000807685986030332, 0.0004731486857077167, 0.00016881480210494433], # 300: [1.431338551855931e-05, 4.012227180120215e-05, 6.850298129947514e-05, 9.605899747347381e-05, 0.00012290072040166636, 0.0001500289025800392, 0.00017760205103018156, 0.0002050376837370814, 0.00023216802273282316, 0.00025936376117528543, 0.000286770923337297, 0.0003141263938546069, 0.0003413064333028032, 0.0003685059563224142, 0.0003958282766766556, 0.00042311960183710317, 0.00045029066797970557, 0.00047746412489335334, 0.0005047156371897132, 0.0005319422306222964, 0.0005590766158027022, 0.0005862043554606536, 0.0006133833418229676, 0.0006405386949363323, 0.0006676176085674061, 0.0006946838416537453, 0.0007217832962185991, 0.0007488584417973781, 0.0007758664924472269, 0.0008028571692832445, 0.0008298678975835217, 0.0008568527124189767, 0.0008837760635199613, 0.0009106780294555469, 0.0009375898032692981, 0.0009644735729342993, 0.000991299201423038, 0.0010180998154806211, 0.001044901876521459, 0.0010716735696885908, 0.0010983889275928184, 0.0011250758736232511, 0.0011517571804946123, 0.00117840559504653, 0.0012049984403109213, 0.0012315596269869425, 0.0012581089880280057, 0.0012846228351818792, 0.0013110811414252545, 0.00133750464674385, 0.001363910797264291, 0.0013902787534588901, 0.0014165906596462637, 0.0014428646992844724, 0.0014691163496578206, 0.0014953270906152946, 0.001521480872306304, 0.0015475937814626148, 0.00157367964899441, 0.0015997218734369169, 0.0016257059263124127, 0.001651646149529894, 0.001677554980842586, 0.001703417427928488, 0.0017292202586412138, 0.0017549763445420545, 0.0017806969321644153, 0.0018063683949693567, 0.0018319786165106394, 0.001857539215647843, 0.0018830604109505785, 0.0019085297473710398, 0.0019339360773297323, 0.001959289942067989, 0.0019846006657902662, 0.002009856807725714, 0.002035048068429775, 0.0020601840541859453, 0.0020852733053521944, 0.00211030526669342, 0.002135270386619155, 0.002160177454021052, 0.0021850343177347107, 0.002209831201507028, 0.0022345592175448837, 0.0022592264352258922, 0.002283840089667751, 0.002308391094568184, 0.0023328711548821417, 0.0023572877027127443, 0.002381647425556331, 0.00240594185202409, 0.002430163219323544, 0.002454318391951272, 0.002478413566344175, 0.0025024408222929406, 0.0025263928773852835, 0.00255027608798504, 0.0025740962081959693, 0.0025978458144661396, 0.0026215180600127255, 0.0026451188441838666, 0.0026686535209823115, 0.0026921151165751785, 0.0027154971809750543, 0.0027388052007325915, 0.0027620441665595044, 0.0027852075136901313, 0.002808289155056237, 0.002831294202880455, 0.0028542273168385085, 0.002877082305828928, 0.002899853416015445, 0.0029225454189358694, 0.002945162671626956, 0.002967699325674426, 0.0029901499343276, 0.0030125189580153013, 0.003034810476250164, 0.0030570189560660717, 0.0030791392346825115, 0.00310117548753865, 0.0031231315389228563, 0.0031450021472680554, 0.0031667824132465515, 0.003188476250471062, 0.003210087247886691, 0.0032316104340022714, 0.0032530411546735326, 0.0032743830823033305, 0.003295639588282603, 0.003316805952230262, 0.003337877748850066, 0.0033588584277646333, 0.003379751158776646, 0.003400551455680898, 0.0034212551073882925, 0.0034418653572642383, 0.0034623851879040167, 0.0034828103321097002, 0.0035031367798385575, 0.0035233675830591412, 0.0035435055501479742, 0.0035635466192892535, 0.0035834869696225067, 0.0036033294751299954, 0.003623076781728268, 0.0036427250207152315, 0.003662270549687044, 0.0036817160767749556, 0.0037010640961099225, 0.0037203109210253195, 0.003739453077859395, 0.0037584931198988143, 0.003777433399208016, 0.003796270401119918, 0.003815000811903878, 0.003833627040006793, 0.003852151304287865, 0.0038705702529859497, 0.003888880724274518, 0.003907084990886804, 0.003925185146571497, 0.003943177994202738, 0.003961060516554395, 0.00397883485897057, 0.00399650299750876, 0.004014061882142001, 0.004031508633574408, 0.0040488452773837025, 0.004066073678737227, 0.004083190927833456, 0.0041001942772090185, 0.004117085639660952, 0.004133866775714363, 0.004150534909511804, 0.004167087419841768, 0.00418352611313258, 0.004199852651012914, 0.004216064385826003, 0.004232158817490263, 0.0042481376519647165, 0.004264002457277361, 0.004279750708703243, 0.004295380022593513, 0.004310892009868351, 0.004326288149837712, 0.004341566035873309, 0.004356723396447817, 0.004371761752443555, 0.004386682498963286, 0.004401483343036293, 0.0044161621212866105, 0.004430720269180384, 0.004445159101775619, 0.004459476435677051, 0.004473670212010881, 0.004487741785091177, 0.0045016923937820415, 0.004515519960510983, 0.0045292225275432845, 0.0045428013719821725, 0.004556257660061908, 0.004569589416572178, 0.004582794781821561, 0.004595874959351039, 0.004608831046059929, 0.004621661165914976, 0.004634363554417776, 0.00464693934491508, 0.00465938956802471, 0.004671712443956493, 0.004683906300772354, 0.0046959722047512755, 0.004707911123050588, 0.0047197213694184096, 0.004731401362055238, 0.004742952103059206, 0.004754374498742801, 0.004765666953897418, 0.004776827974631882, 0.004787858501534513, 0.004798759382490675, 0.004809529111038303, 0.004820166279139213, 0.004830671768345388, 0.004841046370347547, 0.004851288665314363, 0.004861397329170941, 0.004871373186717853, 0.004881216975513267, 0.0048909273604094364, 0.004900503099558432, 0.004909944963118174, 0.0049192536364175265, 0.004928427867200543, 0.004937466494252367, 0.004946370235033431, 0.0049551397243954, 0.004963773791329307, 0.004972271353793974, 0.004980633078343443, 0.004988859550960298, 0.004996949680372135, 0.005004902462380824, 0.005012718514281037, 0.0050203983746564115, 0.005027941030594458, 0.00503534555452279, 0.005042612515983653, 0.005049742407509773, 0.005056734293288795, 0.0050635873212702565, 0.0050703020146247325, 0.005076878821053563, 0.0050833168806995, 0.005089615416040483, 0.005095774905126578, 0.005101795751935285, 0.0051076771715272175, 0.0051134184600114575, 0.005119020051454203, 0.0051244823071031915, 0.005129804516003444, 0.005134986047095221, 0.005140027291483333, 0.005144928568571508, 0.005149689240552303, 0.0051543087484908, 0.005158787441443115, 0.005163125597752833, 0.005167322652015288, 0.0051713781168029205, 0.005175292299931676, 0.005179065439366613, 0.0051826970414525384, 0.005186186689737414, 0.005189534651502356, 0.005192741124918436, 0.005195805687514126, 0.005198727993361184, 0.005201508269816953, 0.0052041466757448375, 0.005206642859375407, 0.005208996544935003, 0.00521120792036838, 0.00521327710562835, 0.005215203819244274, 0.005216987855309183, 0.005218629362770544, 0.005220128422976571, 0.005221484824432338, 0.005222698430884936, 0.005223769352611614, 0.005224697632566832, 0.005225483128989441, 0.005226125775143653, 0.005226625642871415, 0.005226982736856913, 0.005227196984904545, 0.005227268389736952, 0.005227196984904609, 0.005226982736857298, 0.005226625642871094, 0.005226125775143343, 0.005225483128989314, 0.005224697632566667, 0.005223769352611394, 0.005222698430885047, 0.005221484824432276, 0.005220128422976045, 0.005218629362770703, 0.005216987855309369, 0.005215203819244113, 0.005213277105628448, 0.005211207920368179, 0.005208996544935323, 0.005206642859375223, 0.005204146675744738, 0.005201508269816841, 0.005198727993361313, 0.005195805687514172, 0.0051927411249182535, 0.00518953465150263, 0.005186186689737238, 0.0051826970414525254, 0.0051790654393667, 0.005175292299931799, 0.005171378116803026, 0.005167322652015363, 0.005163125597752643, 0.005158787441443411, 0.00515430874849076, 0.005149689240552316, 0.005144928568571719, 0.005140027291483478, 0.00513498604709535, 0.005129804516003328, 0.0051244823071030276, 0.005119020051454297, 0.005113418460011166, 0.005107677171526961, 0.005101795751935227, 0.005095774905126261, 0.005089615416040522, 0.005083316880699735, 0.005076878821053743, 0.005070302014624845, 0.0050635873212700805, 0.0050567342932884085, 0.005049742407509847, 0.00504261251598381, 0.005035345554522823, 0.005027941030594158, 0.005020398374656165, 0.005012718514281176, 0.005004902462380897, 0.00499694968037226, 0.0049888595509603095, 0.004980633078343566, 0.004972271353793644, 0.004963773791329042, 0.004955139724395721, 0.004946370235034016, 0.004937466494252536, 0.004928427867200534, 0.004919253636417707, 0.004909944963118104, 0.004900503099558809, 0.0048909273604097105, 0.0048812169755131024, 0.004871373186717716, 0.004861397329170951, 0.004851288665314311, 0.004841046370347232, 0.004830671768345188, 0.004820166279138959, 0.004809529111038251, 0.004798759382490338, 0.0047878585015345716, 0.004776827974631735, 0.004765666953897307, 0.0047543744987427914, 0.004742952103059202, 0.0047314013620551616, 0.004719721369418516, 0.004707911123050433, 0.004695972204750941, 0.0046839063007723, 0.004671712443956289, 0.004659389568024603, 0.004646939344914928, 0.004634363554417283, 0.00462166116591521, 0.004608831046059741, 0.004595874959351329, 0.004582794781821686, 0.004569589416571941, 0.004556257660061593, 0.004542801371981853, 0.0045292225275430625, 0.0045155199605110725, 0.004501692393782684, 0.004487741785091521, 0.0044736702120110345, 0.0044594764356769125, 0.004445159101775285, 0.0044307202691803455, 0.004416162121286692, 0.004401483343036429, 0.004386682498963625, 0.004371761752443638, 0.004356723396447825, 0.004341566035873333, 0.004326288149837637, 0.004310892009868086, 0.00429538002259381, 0.0042797507087033665, 0.004264002457277777, 0.004248137651964907, 0.004232158817490341, 0.004216064385826206, 0.004199852651012641, 0.004183526113132338, 0.004167087419841528, 0.004150534909512119, 0.004133866775714176, 0.004117085639661458, 0.0041001942772092145, 0.00408319092783341, 0.004066073678737072, 0.004048845277383841, 0.004031508633574328, 0.004014061882142042, 0.003996502997508904, 0.003978834858970341, 0.00396106051655424, 0.003943177994202945, 0.003925185146572065, 0.003907084990886625, 0.0038888807242744516, 0.003870570252985675, 0.0038521513042880833, 0.003833627040007188, 0.0038150008119035845, 0.0037962704011200703, 0.003777433399207321, 0.003758493119898556, 0.003739453077858988, 0.0037203109210250346, 0.0037010640961101984, 0.0036817160767748554, 0.0036622705496868126, 0.00364272502071535, 0.0036230767817284986, 0.003603329475130338, 0.0035834869696225024, 0.0035635466192891576, 0.00354350555014782, 0.0035233675830593134, 0.0035031367798384993, 0.003482810332109164, 0.003462385187904296, 0.003441865357264365, 0.0034212551073889708, 0.003400551455680581, 0.0033797511587767763, 0.003358858427764282, 0.003337877748850245, 0.0033168059522300543, 0.00329563958828269, 0.0032743830823032723, 0.0032530411546732876, 0.003231610434001921, 0.0032100872478863954, 0.003188476250471022, 0.0031667824132466152, 0.0031450021472676755, 0.003123131538923296, 0.003101175487538483, 0.003079139234681824, 0.0030570189560658587, 0.0030348104762501467, 0.003012518958015103, 0.0029901499343269566, 0.0029676993256746486, 0.0029451626716272246, 0.0029225454189359683, 0.002899853416015462, 0.0028770823058289967, 0.0028542273168381985, 0.0028312942028804098, 0.002808289155056215, 0.00278520751368993, 0.0027620441665596666, 0.002738805200732458, 0.002715497180975006, 0.0026921151165756846, 0.0026686535209825175, 0.0026451188441835257, 0.0026215180600122068, 0.002597845814466621, 0.0025740962081962533, 0.0025502760879849754, 0.00252639287738533, 0.002502440822293173, 0.0024784135663443014, 0.002454318391950843, 0.0024301632193235606, 0.0024059418520241804, 0.0023816474255563412, 0.0023572877027127643, 0.002332871154882404, 0.0023083910945681076, 0.002283840089667744, 0.0022592264352261303, 0.0022345592175454193, 0.002209831201507626, 0.0021850343177344657, 0.0021601774540208397, 0.0021352703866189507, 0.0021103052666933977, 0.0020852733053521146, 0.002060184054186014, 0.002035048068429977, 0.002009856807725957, 0.0019846006657902758, 0.00195928994206817, 0.0019339360773296312, 0.0019085297473711007, 0.0018830604109507608, 0.0018575392156482725, 0.0018319786165107877, 0.0018063683949691585, 0.00178069693216468, 0.0017549763445423654, 0.0017292202586409347, 0.0017034174279280823, 0.0016775549808421297, 0.001651646149530077, 0.00162570592631277, 0.0015997218734367048, 0.001573679648994387, 0.0015475937814620998, 0.0015214808723062155, 0.0014953270906152233, 0.0014691163496581476, 0.0014428646992844813, 0.0014165906596466105, 0.0013902787534589855, 0.0013639107972643685, 0.0013375046467437846, 0.0013110811414257816, 0.0012846228351825408, 0.0012581089880279953, 0.0012315596269869729, 0.001204998440311574, 0.001178405595046257, 0.0011517571804948782, 0.0011250758736229627, 0.001098388927592197, 0.001071673569688066, 0.0010449018765211608, 0.0010180998154806183, 0.0009912992014234517, 0.0009644735729341865, 0.0009375898032692254, 0.000910678029455875, 0.0008837760635201092, 0.0008568527124191105, 0.0008298678975834106, 0.0008028571692830543, 0.0007758664924465403, 0.0007488584417971681, 0.000721783296218569, 0.0006946838416533995, 0.0006676176085676541, 0.000640538694936538, 0.0006133833418229547, 0.0005862043554614876, 0.000559076615802472, 0.0005319422306224783, 0.0005047156371895406, 0.0004774641248931814, 0.0004502906679796702, 0.00042311960183722124, 0.00039582827667671103, 0.00036850595632271756, 0.00034130643330288764, 0.0003141263938549482, 0.0002867709233374423, 0.0002593637611755845, 0.0002321680227330003, 0.00020503768373677231, 0.00017760205102980336, 0.00015002890258022111, 0.00012290072040175952, 9.605899747384604e-05, 6.850298129895686e-05, 4.012227180096977e-05, 1.4313385518582305e-05], } kronrod_points = { # 2: [-0.9258200997725493, -0.5773502691896243, 4.440892098500626e-16, 0.5773502691896262, 0.9258200997725514], # 3: [-0.9604912687080189, -0.7745966692414805, -0.4342437493468012, 5.551115123125783e-16, 0.43424374934680354, 0.7745966692414835, 0.9604912687080203], # 4: [-0.9765602507375716, -0.8611363115940519, -0.6402862174963088, -0.3399810435848555, -3.3306690738754696e-16, 0.3399810435848566, 0.6402862174963104, # 0.8611363115940526, 0.9765602507375731 # ], # 5: [-0.9840853600948419, -0.9061798459386637, -0.7541667265708493, -0.5384693101056838, -0.2796304131617827, -5.551115123125783e-16, 0.2796304131617834, # 0.5384693101056831, 0.7541667265708492, 0.9061798459386639, 0.9840853600948425 # ], # 6: [-0.9887032026126787, -0.9324695142031516, -0.821373340865027, -0.6612093864662647, -0.4631182124753035, -0.23861918608319632, - # 2.220446049250313e-16, 0.23861918608319654, 0.4631182124753045, 0.6612093864662645, 0.821373340865028, 0.9324695142031519, 0.9887032026126789 # ], # 7: [-0.9914553711208121, -0.949107912342757, -0.8648644233597671, -0.7415311855993938, -0.5860872354676903, -0.40584515137739663, -0.20778495500789784, # 3.3306690738754696e-16, 0.20778495500789873, 0.40584515137739796, 0.5860872354676913, 0.7415311855993948, 0.8648644233597691, # 0.9491079123427586, 0.9914553711208126 # ], # 8: [-0.9933798758817153, -0.9602898564975361, -0.8941209068474557, -0.7966664774136257, -0.6723540709451575, -0.5255324099163289, -0.3607010979281312, - # 0.18343464249564911, 4.440892098500626e-16, 0.1834346424956501, 0.3607010979281319, 0.5255324099163294, 0.6723540709451585, 0.7966664774136267, # 0.8941209068474564, 0.9602898564975363, 0.9933798758817162 # ], # 9: [-0.9946781606773386, -0.9681602395076251, -0.9149635072496775, -0.8360311073266354, -0.7344867651839325, -0.6133714327005891, -0.475462479112459, - # 0.3242534234038089, -0.16422356361498636, 3.3306690738754696e-16, 0.16422356361498702, 0.3242534234038096, 0.4754624791124603, # 0.613371432700591, 0.7344867651839341, 0.8360311073266359, 0.914963507249678, 0.9681602395076261, 0.9946781606773403 # ], 10: [-0.9956571630258076, -0.9739065285171714, -0.9301574913557086, -0.8650633666889832, -0.7808177265864168, -0.679409568299023, -0.562757134668605, - 0.4333953941292482, -0.29439286270146015, -0.14887433898163127, 2.220446049250313e-16, 0.1488743389816316, 0.2943928627014607, 0.43339539412924755, 0.5627571346686051, 0.6794095682990244, 0.7808177265864169, 0.8650633666889844, 0.9301574913557082, 0.9739065285171717, 0.9956571630258081 ], # 11: [-0.9963696138895425, -0.9782286581460574, -0.9416771085780676, -0.8870625997680958, -0.8160574566562212, -0.7301520055740492, -0.6305995201619655, - # 0.5190961292068118, -0.39794414095237773, -0.26954315595234524, -0.13611300079936217, -4.440892098500626e-16, 0.13611300079936206, # 0.2695431559523449, 0.39794414095237773, 0.5190961292068117, 0.6305995201619654, 0.7301520055740492, 0.816057456656221, 0.8870625997680953, # 0.9416771085780679, 0.978228658146057, 0.9963696138895426 # ], # 12: [-0.9969339225295953, -0.9815606342467201, -0.950537795943121, -0.9041172563704741, -0.8435581241611524, -0.7699026741943046, -0.6840598954700551, - # 0.5873179542866177, -0.4813394504781564, -0.3678314989981798, -0.24850574832046946, -0.12523340851146914, -3.3306690738754696e-16, # 0.1252334085114689, 0.24850574832046957, 0.36783149899818013, 0.4813394504781572, 0.5873179542866174, 0.6840598954700562, 0.7699026741943048, # 0.8435581241611532, 0.9041172563704747, 0.9505377959431213, 0.9815606342467192, 0.9969339225295955 # ], # 13: [-0.9973661769948244, -0.984183054718588, -0.95755246838608, -0.9175983992229771, -0.8653331602663435, -0.8015780907333098, -0.7269488493206314, - # 0.64234933944034, -0.5490799579565359, -0.4484927510364468, -0.34183246302180625, -0.23045831595513466, -0.11597108974493386, - # 3.3306690738754696e-16, 0.11597108974493409, 0.2304583159551351, 0.3418324630218067, 0.44849275103644715, 0.5490799579565372, # 0.6423493394403403, 0.726948849320632, 0.80157809073331, 0.8653331602663444, 0.917598399222978, 0.9575524683860812, 0.9841830547185881, # 0.9973661769948249 # ], # 14: [-0.9977205937565432, -0.9862838086968115, -0.9631583382788527, -0.9284348836635736, -0.8829146632520561, -0.8272013150697641, -0.7617567525622055, - # 0.6872929048116849, -0.6047893659409216, -0.5152486363581534, -0.41965589764297806, -0.31911236892788875, -0.21483591853348383, - # 0.108054948707343, 5.551115123125783e-16, 0.10805494870734478, 0.21483591853348527, 0.3191123689278901, 0.41965589764297917, 0.5152486363581541, # 0.6047893659409218, 0.6872929048116857, 0.7617567525622057, 0.8272013150697651, 0.8829146632520571, 0.9284348836635736, 0.9631583382788532, # 0.9862838086968123, 0.9977205937565431 # ], # 15: [-0.9980022986933965, -0.9879925180204847, -0.9677390756791384, -0.937273392400706, -0.8972645323440815, -0.8482065834104267, -0.7904185014424661, - # 0.7244177313601698, -0.6509967412974171, -0.5709721726085399, -0.48508186364023964, -0.3941513470775635, -0.2991800071531687, - # 0.20119409399743393, -0.1011420669187173, 1.1102230246251565e-16, 0.10114206691871763, 0.20119409399743426, 0.2991800071531686, # 0.3941513470775634, 0.48508186364023986, 0.570972172608539, 0.6509967412974169, 0.7244177313601701, 0.7904185014424661, 0.8482065834104272, # 0.897264532344082, 0.937273392400706, 0.9677390756791392, 0.9879925180204855, 0.9980022986933971 # ], # 16: [-0.9982392741454447, -0.9894009349916497, -0.9715059509693929, -0.9445750230732326, -0.9091576670123429, -0.8656312023878314, -0.8142402870624446, - # 0.7554044083550029, -0.6897411066817623, -0.6178762444026438, -0.5404076763521388, -0.45801677765722726, -0.37148378087841594, - # 0.28160355077925914, -0.18916857901808393, -0.09501250983763732, -2.220446049250313e-16, 0.0950125098376372, 0.18916857901808315, # 0.2816035507792588, 0.3714837808784164, 0.45801677765722737, 0.5404076763521395, 0.6178762444026438, 0.6897411066817624, 0.7554044083550031, # 0.8142402870624446, 0.8656312023878316, 0.9091576670123429, 0.9445750230732325, 0.9715059509693925, 0.9894009349916499, 0.9982392741454444 # ], # 17: [-0.998432970606058, -0.9905754753144171, -0.9746592569674308, -0.9506755217687672, -0.9190961368038919, -0.8802391537269854, -0.8342740928501335, - # 0.7815140038968009, -0.722472287372409, -0.6576711592166902, -0.587569212334035, -0.5126905370864765, -0.43368729520979854, -0.3512317634538762, - # 0.2659465074516818, -0.17848418149584755, -0.08958563942522657, 0.0, 0.0895856394252269, 0.17848418149584822, 0.26594650745168247, # 0.3512317634538764, 0.43368729520979965, 0.5126905370864772, 0.5875692123340355, 0.6576711592166908, 0.7224722873724099, 0.7815140038968015, # 0.8342740928501343, 0.8802391537269859, 0.9190961368038917, 0.9506755217687678, 0.974659256967431, 0.9905754753144174, 0.9984329706060581 # ], # 18: [-0.998599165412682, -0.9915651684209306, -0.9773103807748602, -0.9558239495713969, -0.927504501573343, -0.8926024664975549, -0.8512493734337165, - # 0.8037049589725229, -0.7503804817033572, -0.6916870430603526, -0.6280037548638335, -0.5597708310739471, -0.48750904258507766, - # 0.4117511614628424, -0.3330234227914082, -0.2518862256915051, -0.16893682806654864, -0.08477501304173507, 0.0, 0.0847750130417354, # 0.1689368280665493, 0.2518862256915053, 0.33302342279140795, 0.41175116146284263, 0.48750904258507877, 0.5597708310739478, 0.6280037548638338, # 0.6916870430603532, 0.7503804817033577, 0.803704958972523, 0.8512493734337176, 0.8926024664975558, 0.9275045015733433, 0.9558239495713977, # 0.977310380774861, 0.9915651684209309, 0.9985991654126828 # ], # 19: [-0.9987380120803147, -0.9924068438435842, -0.9795723007892643, -0.9602081521348296, -0.9346640001655193, -0.903155903614818, -0.865773542685335, - # 0.822714656537143, -0.7743288590095345, -0.7209661773352292, -0.6629230657449079, -0.6005453046616815, -0.5342756325305535, - # 0.46457074137596144, -0.3918512308715205, -0.3165640999636299, -0.23922492831824171, -0.16035864564022484, -0.0804519227440117, - # 3.3306690738754696e-16, 0.08045192274401136, 0.16035864564022517, 0.23922492831824105, 0.31656409996363, 0.3918512308715202, 0.464570741375961, # 0.5342756325305533, 0.6005453046616811, 0.6629230657449078, 0.7209661773352294, 0.7743288590095347, 0.8227146565371429, 0.8657735426853352, # 0.9031559036148178, 0.9346640001655193, 0.9602081521348299, 0.9795723007892634, 0.9924068438435844, 0.9987380120803145 # ], # 87: [-0.9999373393176647, -0.999622350432484, -0.998981834567236, -0.9980107215629824, -0.9967203521912877, -0.9951134581379919, -0.9931850143689511, -0.9909336501109647, -0.9883633935160456, -0.985476606068421, -0.9822719826624027, -0.978749341003817, -0.9749114465404709, -0.9707605197472842, -0.9662966427216388, -0.961520436936073, -0.9564344462027763, -0.9510410013893145, -0.945341012434413, -0.9393357199742063, -0.9330277577952324, -0.9264196798912376, -0.9195130155009865, -0.9123095290941232, -0.9048120659287133, -0.8970234547675697, -0.8889457425315686, -0.8805811598520423, -0.8719328161808033, -0.8630038376297071, -0.8537967282834917, -0.8443141443970357, -0.8345594815785893, -0.8245361702557805, -0.8142471292361737, -0.8036954080581786, -0.7928846943811525, -0.7818187205454386, -0.7705007861856665, -0.7589343057213496, -0.7471232552715923, -0.735071660505314, -0.7227831721000328, -0.7102615427114733, -0.6975110271436125, -0.6845359314028241, -0.6713402293006723, -0.6579279856713498, -0.6443037198655162, -0.6304720018972505, -0.6164371014898609, -0.6022033695423357, -0.5877755725413023, -0.5731585255345304, -0.5583567669734248, -0.5433749073022942, -0.5282179401678413, -0.5128909045194527, -0.49739858004223114, -0.4817458096226672, -0.4659377920214407, -0.449979767165285, -0.43387672799250576, -0.417633722072869, -0.4012561295371335, -0.3847493668670974, -0.3681186117187589, -0.3513690879260418, -0.33450633184746614, -0.3175359108501159, -0.30046315821599157, -0.28329344500409226, -0.266032437514764, -0.24868582730667765, -0.2312590736832678, -0.21375766534010188, -0.19618737131686093, -0.17855397985304466, -0.16086304622927927, -0.14312014673100815, -0.12533112522788015, -0.10750183850894302, -0.08963790744155298, -0.07174496550149989, -0.05382890297085785, -0.03589561662680141, -0.017950762293751588, 0.0, 0.017950762293751477, 0.03589561662680074, 0.05382890297085752, 0.071744965501499, 0.08963790744155298, 0.10750183850894268, 0.1253311252278798, 0.14312014673100815, 0.16086304622927938, 0.17855397985304444, 0.19618737131686126, 0.21375766534010154, 0.23125907368326792, 0.2486858273066771, 0.2660324375147639, 0.28329344500409204, 0.30046315821599134, 0.3175359108501157, 0.3345063318474658, 0.3513690879260414, 0.3681186117187585, 0.38474936686709693, 0.4012561295371332, 0.4176337220728691, 0.4338767279925052, 0.4499797671652843, 0.46593779202144003, 0.4817458096226668, 0.49739858004223125, 0.5128909045194523, 0.5282179401678414, 0.5433749073022933, 0.5583567669734248, 0.5731585255345304, 0.5877755725413025, 0.6022033695423353, 0.6164371014898609, 0.6304720018972505, 0.6443037198655157, 0.6579279856713494, 0.671340229300672, 0.6845359314028239, 0.6975110271436127, 0.7102615427114729, 0.7227831721000324, 0.7350716605053138, 0.7471232552715918, 0.7589343057213493, 0.7705007861856658, 0.7818187205454378, 0.7928846943811524, 0.8036954080581777, 0.8142471292361735, 0.8245361702557805, 0.8345594815785895, 0.8443141443970352, 0.8537967282834915, 0.8630038376297067, 0.8719328161808031, 0.8805811598520419, 0.8889457425315681, 0.8970234547675693, 0.9048120659287131, 0.9123095290941229, 0.9195130155009866, 0.9264196798912372, 0.933027757795232, 0.9393357199742057, 0.9453410124344126, 0.9510410013893136, 0.9564344462027765, 0.9615204369360731, 0.9662966427216388, 0.9707605197472843, 0.9749114465404707, 0.9787493410038163, 0.9822719826624021, 0.9854766060684204, 0.988363393516045, 0.9909336501109646, 0.9931850143689511, 0.995113458137992, 0.9967203521912877, 0.9980107215629823, 0.998981834567236, 0.9996223504324843, 0.9999373393176647], # 300: [-0.9999946872996457, -0.9999679782184366, -0.9999136584901585, -0.9998312829844194, -0.9997217884834974, -0.9995853734488863, -0.9994215723741114, -0.9992302218632736, -0.9990116062590783, -0.9987658603530949, -0.9984928023779717, -0.9981923380734435, -0.9978646126899144, -0.9975097171758793, -0.9971275560703853, -0.996718072051, -0.9962813600745217, -0.9958174891208241, -0.9953264031603122, -0.9948080667630033, -0.9942625517027517, -0.9936899152742682, -0.9930901239339065, -0.9924631568471005, -0.9918090737598455, -0.9911279255509896, -0.990419693459469, -0.9896843673150096, -0.9889220003217568, -0.9881326399103605, -0.9873162780530513, -0.9864729129322606, -0.9856025942467834, -0.9847053677808956, -0.9837812338860811, -0.9828301976412723, -0.9818523070673794, -0.9808476074619001, -0.9798161060814252, -0.9787578139322586, -0.9776727785275657, -0.9765610454590224, -0.9754226279104833, -0.9742575421410374, -0.9730658359415845, -0.9718475557430342, -0.9706027199771077, -0.9693313496680812, -0.9680334934315086, -0.9667091989287729, -0.9653584893488589, -0.9639813901171553, -0.9625779510657028, -0.9611482213733913, -0.9596922286202931, -0.9582100023531643, -0.9567015939075111, -0.9551670541938286, -0.9536064149019075, -0.9520197094793031, -0.9504069909786793, -0.9487683122037178, -0.9471037087320028, -0.9454132177338136, -0.9436968941399522, -0.9419547927700981, -0.9401869529103559, -0.9383934153067652, -0.9365742368903639, -0.9347294745903799, -0.9328591712533739, -0.9309633710773335, -0.929042133086303, -0.9270955163900743, -0.9251235672708156, -0.9231263332721165, -0.9211038755812251, -0.9190562555418369, -0.9169835227643917, -0.9148857280450582, -0.9127629347867887, -0.9106152066064295, -0.9084425963487075, -0.9062451579796027, -0.9040229571561673, -0.901776059796223, -0.8995045218950815, -0.8972084005137191, -0.894887763590253, -0.8925426793619207, -0.8901732068988839, -0.8877794062884896, -0.8853613477675033, -0.8829191019031915, -0.8804527307710444, -0.8779622974209855, -0.8754478743981893, -0.8729095346039717, -0.8703473430544713, -0.8677613657021737, -0.8651516774038297, -0.8625183533931811, -0.8598614615661291, -0.8571810707206511, -0.8544772580226088, -0.8517501010316768, -0.848999670465588, -0.846226037913026, -0.8434292828416308, -0.8406094851262689, -0.8377667182508701, -0.8349010565417854, -0.832012581756866, -0.8291013760716741, -0.8261675156824684, -0.8232110776015575, -0.8202321458616895, -0.8172308049213011, -0.8142071336364427, -0.8111612116546523, -0.8080931252649417, -0.8050029611877987, -0.8018908008875274, -0.7987567265968607, -0.7956008268394583, -0.7924231905743322, -0.7892239018232119, -0.786003045354458, -0.7827607119020659, -0.7794969926375797, -0.7762119740898052, -0.7729057435134399, -0.7695783938260492, -0.7662300183834612, -0.7628607061714925, -0.7594705468820133, -0.7560596355871452, -0.752628067796661, -0.7491759349034655, -0.7457033289874226, -0.7422103472441275, -0.7386970873050263, -0.7351636429201902, -0.7316101085081914, -0.7280365833550897, -0.7244431671799338, -0.7208299560399418, -0.7171970466429054, -0.7135445403305545, -0.709872538873781, -0.706181140586744, -0.7024704444166954, -0.6987405537245666, -0.69499157229575, -0.6912236006508845, -0.6874367399265875, -0.6836310954649563, -0.679806773027049, -0.675963875289201, -0.6721025055269298, -0.6682227710239814, -0.6643247794768312, -0.6604086356663699, -0.6564744449561277, -0.6525223165305882, -0.6485523599800562, -0.6445646821384352, -0.640559390405937, -0.6365365958255456, -0.632496409838541, -0.6284389412798634, -0.6243642995346044, -0.6202725974607572, -0.6161639483065122, -0.612038462855417, -0.6078962524251541, -0.6037374316440418, -0.5995621155219637, -0.5953704167381737, -0.591162448490157, -0.5869383271307362, -0.5826981693851709, -0.5784420897750377, -0.5741702033243279, -0.569882628063465, -0.5655794823857194, -0.5612608826011143, -0.5569269455063359, -0.552577790761474, -0.5482135383794409, -0.5438343064043345, -0.539440213351214, -0.5350313804609139, -0.5306079293166595, -0.5261699796417558, -0.5217176516148028, -0.5172510680075144, -0.512770351923185, -0.5082756247089657, -0.5037670081516556, -0.4992446265030859, -0.4947086043354987, -0.49015906456404057, -0.48559613052787265, -0.4810199279073166, -0.47643058269160643, -0.4718282193075505, -0.4672129625903425, -0.4625849395963533, -0.4579442776790468, -0.45329110272009054, -0.4486255409938984, -0.4439477208796667, -0.43925777104156594, -0.4345558187588534, -0.4298419916878913, -0.42511641947672607, -0.42037923204598227, -0.415630558014791, -0.4108705263637342, -0.4060992679550248, -0.40131691391079194, -0.3965235941318921, -0.39171943886495786, -0.38690458013101003, -0.3820791501980705, -0.3772432801901535, -0.37239710156135253, -0.36754074743549325, -0.36267435117025193, -0.3577980450538263, -0.35291196168878125, -0.3480162352451315, -0.34311100011337237, -0.338196389686529, -0.3332725376562593, -0.3283395791815773, -0.323397649628395, -0.31844688343484195, -0.31348741532192625, -0.30851938137992496, -0.3035429178922302, -0.29855816028200666, -0.29356524423952535, -0.28856430672807254, -0.28355548489008364, -0.27853891507336437, -0.27351473387703384, -0.26848307907865343, -0.2634440886207914, -0.258397899715193, -0.25334464980911076, -0.24828447743519133, -0.2432175212767862, -0.23814391934859436, -0.23306380988500508, -0.22797733211414206, -0.22288462540037246, -0.21778582850010708, -0.21268108037361766, -0.20757052088450112, -0.2024542900179942, -0.1973325272107358, -0.19220537208739352, -0.18707296508667115, -0.1819354467541754, -0.17679295714507726, -0.17164563648674025, -0.16649362573227822, -0.1613370659268425, -0.15617609768225016, -0.15101086176667544, -0.14584149958664194, -0.1406681526257254, -0.13549096199033794, -0.13031006892741814, -0.12512561523562193, -0.1199377427755659, -0.11474659308606205, -0.10955230783063419, -0.10435502913854622, -0.09915489918583553, -0.09395205988140298, -0.08874665324306696, -0.0835388216689601, -0.07832870758870736, -0.07311645321890214, -0.06790220086927956, -0.06268609314491269, -0.05746827266700372, -0.0522488818973752, -0.04702806337524312, -0.04180595985053315, -0.03658271407386193, -0.03135846868977743, -0.026133366404510694, -0.020907550050621282, -0.01568116244584883, -0.010454346354956456, -0.0052272445887175945, 0.0, 0.0052272445887177055, 0.0104543463549569, 0.01568116244584905, 0.020907550050621393, 0.026133366404511027, 0.03135846868977776, 0.036582714073862044, 0.04180595985053359, 0.04702806337524368, 0.052248881897375865, 0.057468272667004605, 0.06268609314491314, 0.0679022008692799, 0.07311645321890248, 0.07832870758870769, 0.08353882166896032, 0.08874665324306752, 0.09395205988140376, 0.09915489918583587, 0.10435502913854666, 0.10955230783063508, 0.1147465930860625, 0.11993774277556657, 0.12512561523562216, 0.13031006892741892, 0.13549096199033805, 0.14066815262572618, 0.14584149958664216, 0.1510108617666761, 0.15617609768225016, 0.16133706592684272, 0.16649362573227855, 0.17164563648674003, 0.17679295714507792, 0.18193544675417572, 0.18707296508667137, 0.1922053720873934, 0.19733252721073624, 0.2024542900179942, 0.20757052088450112, 0.21268108037361744, 0.21778582850010697, 0.22288462540037235, 0.22797733211414206, 0.23306380988500508, 0.23814391934859414, 0.2432175212767861, 0.24828447743519155, 0.25334464980911076, 0.2583978997151932, 0.2634440886207915, 0.2684830790786531, 0.27351473387703407, 0.27853891507336426, 0.2835554848900833, 0.2885643067280719, 0.29356524423952524, 0.2985581602820062, 0.30354291789222965, 0.30851938137992485, 0.3134874153219265, 0.31844688343484184, 0.323397649628395, 0.3283395791815772, 0.33327253765625897, 0.33819638968652843, 0.3431110001133719, 0.34801623524513126, 0.3529119616887808, 0.3577980450538262, 0.3626743511702518, 0.36754074743549336, 0.37239710156135264, 0.3772432801901534, 0.38207915019807037, 0.3869045801310097, 0.39171943886495786, 0.3965235941318921, 0.4013169139107917, 0.4060992679550246, 0.410870526363734, 0.4156305580147909, 0.4203792320459817, 0.42511641947672574, 0.42984199168789106, 0.43455581875885274, 0.4392577710415656, 0.44394772087966616, 0.44862554099389795, 0.45329110272008966, 0.45794427767904666, 0.46258493959635283, 0.4672129625903423, 0.4718282193075505, 0.4764305826916061, 0.4810199279073162, 0.4855961305278721, 0.49015906456404, 0.4947086043354981, 0.49924462650308543, 0.5037670081516556, 0.508275624708965, 0.5127703519231843, 0.517251068007514, 0.5217176516148028, 0.5261699796417558, 0.5306079293166588, 0.5350313804609135, 0.5394402133512137, 0.5438343064043343, 0.5482135383794404, 0.5525777907614735, 0.5569269455063353, 0.5612608826011143, 0.5655794823857192, 0.5698826280634648, 0.5741702033243274, 0.578442089775038, 0.5826981693851705, 0.5869383271307357, 0.5911624484901564, 0.5953704167381735, 0.5995621155219635, 0.6037374316440416, 0.6078962524251535, 0.6120384628554165, 0.6161639483065119, 0.6202725974607567, 0.6243642995346039, 0.6284389412798628, 0.6324964098385406, 0.6365365958255453, 0.6405593904059368, 0.6445646821384348, 0.6485523599800556, 0.6525223165305876, 0.6564744449561274, 0.6604086356663694, 0.6643247794768309, 0.6682227710239813, 0.6721025055269294, 0.6759638752892005, 0.6798067730270486, 0.6836310954649557, 0.687436739926587, 0.6912236006508838, 0.69499157229575, 0.6987405537245661, 0.7024704444166952, 0.706181140586744, 0.7098725388737805, 0.713544540330554, 0.7171970466429054, 0.7208299560399416, 0.7244431671799337, 0.7280365833550895, 0.7316101085081913, 0.7351636429201899, 0.7386970873050263, 0.742210347244127, 0.7457033289874225, 0.7491759349034652, 0.752628067796661, 0.7560596355871447, 0.7594705468820128, 0.7628607061714925, 0.7662300183834609, 0.7695783938260488, 0.7729057435134394, 0.7762119740898051, 0.77949699263758, 0.7827607119020656, 0.7860030453544584, 0.7892239018232124, 0.7924231905743324, 0.7956008268394587, 0.7987567265968609, 0.8018908008875274, 0.8050029611877982, 0.8080931252649417, 0.8111612116546523, 0.8142071336364427, 0.8172308049213008, 0.8202321458616895, 0.8232110776015571, 0.826167515682468, 0.8291013760716742, 0.8320125817568655, 0.8349010565417858, 0.8377667182508699, 0.8406094851262688, 0.8434292828416307, 0.846226037913026, 0.848999670465588, 0.8517501010316768, 0.8544772580226087, 0.8571810707206517, 0.859861461566129, 0.8625183533931813, 0.8651516774038295, 0.8677613657021737, 0.870347343054471, 0.8729095346039718, 0.8754478743981893, 0.8779622974209857, 0.8804527307710446, 0.8829191019031915, 0.8853613477675031, 0.88777940628849, 0.890173206898884, 0.8925426793619202, 0.8948877635902533, 0.8972084005137192, 0.8995045218950817, 0.9017760597962232, 0.9040229571561679, 0.9062451579796031, 0.9084425963487072, 0.9106152066064294, 0.9127629347867888, 0.9148857280450583, 0.9169835227643919, 0.9190562555418367, 0.9211038755812248, 0.9231263332721162, 0.9251235672708154, 0.927095516390074, 0.9290421330863031, 0.9309633710773336, 0.932859171253374, 0.93472947459038, 0.9365742368903638, 0.9383934153067653, 0.940186952910356, 0.9419547927700982, 0.9436968941399521, 0.9454132177338137, 0.9471037087320029, 0.948768312203718, 0.9504069909786793, 0.9520197094793031, 0.9536064149019075, 0.9551670541938287, 0.9567015939075112, 0.9582100023531646, 0.9596922286202932, 0.9611482213733913, 0.9625779510657029, 0.9639813901171553, 0.965358489348859, 0.9667091989287728, 0.9680334934315088, 0.9693313496680812, 0.9706027199771078, 0.9718475557430342, 0.9730658359415845, 0.9742575421410374, 0.9754226279104833, 0.9765610454590226, 0.9776727785275656, 0.9787578139322587, 0.9798161060814253, 0.9808476074619004, 0.9818523070673794, 0.9828301976412724, 0.9837812338860812, 0.9847053677808958, 0.9856025942467834, 0.9864729129322606, 0.9873162780530514, 0.9881326399103604, 0.9889220003217569, 0.9896843673150097, 0.9904196934594691, 0.9911279255509898, 0.9918090737598456, 0.9924631568471006, 0.9930901239339065, 0.9936899152742682, 0.9942625517027519, 0.9948080667630034, 0.9953264031603122, 0.9958174891208242, 0.9962813600745218, 0.9967180720510002, 0.9971275560703854, 0.9975097171758793, 0.9978646126899146, 0.9981923380734435, 0.9984928023779718, 0.998765860353095, 0.9990116062590781, 0.9992302218632736, 0.9994215723741117, 0.9995853734488863, 0.9997217884834975, 0.9998312829844194, 0.9999136584901583, 0.9999679782184367, 0.9999946872996458], } legendre_weights = { # 2: [1.0, 1.0], # 3: [0.5555555555555557, 0.8888888888888888, 0.5555555555555557], # 4: [0.3478548451374537, 0.6521451548625462, 0.6521451548625462, 0.3478548451374537], # 5: [0.23692688505618942, 0.4786286704993662, 0.568888888888889, 0.4786286704993662, 0.23692688505618942], # 6: [0.17132449237916975, 0.36076157304813894, 0.46791393457269137, 0.46791393457269137, 0.36076157304813894, 0.17132449237916975], # 7: [0.12948496616887065, 0.2797053914892766, 0.3818300505051183, 0.41795918367346896, 0.3818300505051183, 0.2797053914892766, 0.12948496616887065], # 8: [0.10122853629037669, 0.22238103445337434, 0.31370664587788705, 0.36268378337836177, 0.36268378337836177, 0.31370664587788705, 0.22238103445337434, # 0.10122853629037669 # ], # 9: [0.08127438836157472, 0.18064816069485712, 0.26061069640293566, 0.3123470770400028, 0.33023935500125967, 0.3123470770400028, 0.26061069640293566, # 0.18064816069485712, 0.08127438836157472 # ], 10: [0.06667134430868807, 0.14945134915058036, 0.219086362515982, 0.2692667193099965, 0.295524224714753, 0.295524224714753, 0.2692667193099965, 0.219086362515982, 0.14945134915058036, 0.06667134430868807 ], # 11: [0.055668567116173164, 0.1255803694649047, 0.18629021092773443, 0.23319376459199068, 0.26280454451024676, 0.2729250867779009, 0.26280454451024676, # 0.23319376459199068, 0.18629021092773443, 0.1255803694649047, 0.055668567116173164 # ], # 12: [0.04717533638651202, 0.10693932599531888, 0.1600783285433461, 0.20316742672306565, 0.23349253653835464, 0.2491470458134027, 0.2491470458134027, # 0.23349253653835464, 0.20316742672306565, 0.1600783285433461, 0.10693932599531888, 0.04717533638651202 # ], # 13: [0.04048400476531588, 0.0921214998377286, 0.13887351021978736, 0.17814598076194552, 0.20781604753688857, 0.22628318026289715, 0.2325515532308739, # 0.22628318026289715, 0.20781604753688857, 0.17814598076194552, 0.13887351021978736, 0.0921214998377286, 0.04048400476531588 # ], # 14: [0.035119460331752374, 0.0801580871597603, 0.12151857068790296, 0.1572031671581934, 0.18553839747793763, 0.20519846372129555, 0.21526385346315766, # 0.21526385346315766, 0.20519846372129555, 0.18553839747793763, 0.1572031671581934, 0.12151857068790296, 0.0801580871597603, # 0.035119460331752374 # ], # 15: [0.030753241996118647, 0.07036604748810807, 0.10715922046717177, 0.1395706779261539, 0.16626920581699378, 0.18616100001556188, 0.19843148532711125, # 0.2025782419255609, 0.19843148532711125, 0.18616100001556188, 0.16626920581699378, 0.1395706779261539, 0.10715922046717177, 0.07036604748810807, # 0.030753241996118647 # ], # 16: [0.027152459411754037, 0.062253523938647706, 0.09515851168249259, 0.12462897125553403, 0.14959598881657676, 0.16915651939500262, 0.1826034150449236, # 0.18945061045506859, 0.18945061045506859, 0.1826034150449236, 0.16915651939500262, 0.14959598881657676, 0.12462897125553403, # 0.09515851168249259, 0.062253523938647706, 0.027152459411754037 # ], # 17: [0.02414830286854952, 0.0554595293739866, 0.08503614831717908, 0.11188384719340365, 0.13513636846852523, 0.15404576107681012, 0.16800410215644995, # 0.17656270536699253, 0.17944647035620653, 0.17656270536699253, 0.16800410215644995, 0.15404576107681012, 0.13513636846852523, # 0.11188384719340365, 0.08503614831717908, 0.0554595293739866, 0.02414830286854952 # ], # 18: [0.02161601352648413, 0.04971454889496922, 0.07642573025488925, 0.10094204410628699, 0.12255520671147836, 0.14064291467065063, 0.15468467512626521, # 0.16427648374583273, 0.16914238296314363, 0.16914238296314363, 0.16427648374583273, 0.15468467512626521, 0.14064291467065063, # 0.12255520671147836, 0.10094204410628699, 0.07642573025488925, 0.04971454889496922, 0.02161601352648413 # ], # 19: [0.01946178822972761, 0.04481422676569981, 0.06904454273764107, 0.09149002162244985, 0.11156664554733375, 0.1287539625393362, 0.14260670217360638, # 0.15276604206585945, 0.15896884339395415, 0.16105444984878345, 0.15896884339395415, 0.15276604206585945, 0.14260670217360638, # 0.1287539625393362, 0.11156664554733375, 0.09149002162244985, 0.06904454273764107, 0.04481422676569981, 0.01946178822972761 # ], # 87: [0.0009691097381757348, 0.0022546907537549384, 0.003539271655388117, 0.004819456238501684, 0.00609346204763528, 0.007359623648817653, 0.008616302838488783, 0.009861877713702014, 0.011094741940560695, 0.012313306030048123, 0.013515999118245947, 0.0147012708872398, 0.015867593518826113, 0.017013463643001287, 0.0181374042653544, 0.01923796666535654, 0.020313732260655384, 0.021363314433802724, 0.022385360318485367, 0.02337855254266003, 0.024341610926167583, 0.025273294130557043, 0.026172401258933563, 0.027037773403735973, 0.02786829514042919, 0.028662895965176235, 0.02942055167462312, 0.03014028568601895, 0.030821170295962166, 0.031462327876150824, 0.032062932004589706, 0.03262220853080143, 0.03313943657366203, 0.03361394945057692, 0.034045135536799345, 0.03443243905378225, 0.03477536078554787, 0.03507345872215154, 0.03532634862941022, 0.03553370454416061, 0.035695259194409426, 0.03581080434383377, 0.03588019106018702, 0.03590332990726553, 0.03588019106018702, 0.03581080434383377, 0.035695259194409426, 0.03553370454416061, 0.03532634862941022, 0.03507345872215154, 0.03477536078554787, 0.03443243905378225, 0.034045135536799345, 0.03361394945057692, 0.03313943657366203, 0.03262220853080143, 0.032062932004589706, 0.031462327876150824, 0.030821170295962166, 0.03014028568601895, 0.02942055167462312, 0.028662895965176235, 0.02786829514042919, 0.027037773403735973, 0.026172401258933563, 0.025273294130557043, 0.024341610926167583, 0.02337855254266003, 0.022385360318485367, 0.021363314433802724, 0.020313732260655384, 0.01923796666535654, 0.0181374042653544, 0.017013463643001287, 0.015867593518826113, 0.0147012708872398, 0.013515999118245947, 0.012313306030048123, 0.011094741940560695, 0.009861877713702014, 0.008616302838488783, 0.007359623648817653, 0.00609346204763528, 0.004819456238501684, 0.003539271655388117, 0.0022546907537549384, 0.0009691097381757348], # 300: [8.217779368318766e-05, 0.000191285544654933, 0.00030053396541231163, 0.0004097635734033533, 0.0005189512134602094, 0.0006280829779724268, 0.0007371464159049574, 0.0008461294290789873, 0.0009550200344351969, 0.0010638062980465963, 0.001172476313703897, 0.001281018195428937, 0.0013894200749916825, 0.0014976701014677357, 0.0016057564416445678, 0.0017136672808561744, 0.0018213908240185923, 0.0019289152967680268, 0.0020362289466670158, 0.0021433200444213973, 0.002250176885138956, 0.002356787789583519, 0.002463141105430975, 0.0025692252085409055, 0.002675028504212449, 0.002780539428452415, 0.0028857464492311412, 0.0029906380677445906, 0.0030952028196669814, 0.0031994292764019384, 0.003303306046330757, 0.003406821776058694, 0.0035099651516518333, 0.0036127248998770794, 0.0037150897894316927, 0.0038170486321696183, 0.003918590284326712, 0.004019703647735881, 0.004120377671042578, 0.004220601350909488, 0.004320363733222272, 0.0044196539142855205, 0.004518461042012614, 0.004616774317114607, 0.004714582994277626, 0.004811876383340857, 0.004908643850461695, 0.005004874819279474, 0.0051005587720714925, 0.0051956852509013034, 0.005290243858763116, 0.0053842242607182395, 0.0054776161850220466, 0.005570409424251355, 0.005662593836414576, 0.005754159346065033, 0.005845095945399296, 0.005935393695350927, 0.006025042726679271, 0.006114033241045046, 0.006202355512083302, 0.0062899998864659286, 0.00637695678495604, 0.006463216703456925, 0.006548770214047811, 0.006633607966017162, 0.006717720686882713, 0.0068010991834066194, 0.006883734342598478, 0.0069656171327117265, 0.00704673860423293, 0.007127089890856632, 0.007206662210456636, 0.0072854468660452166, 0.007363435246723478, 0.00744061882862333, 0.007516989175837821, 0.00759253794134379, 0.007667256867914979, 0.007741137789022482, 0.007814172629730136, 0.007886353407574561, 0.007957672233438363, 0.008028121312413226, 0.008097692944651201, 0.008166379526205846, 0.008234173549864183, 0.008301067605967079, 0.008367054383218206, 0.008432126669483946, 0.00849627735258214, 0.00855949942105765, 0.0086217859649506, 0.00868313017655067, 0.008743525351140993, 0.008802964887731923, 0.008861442289781127, 0.008918951165904783, 0.008975485230575371, 0.009031038304809223, 0.009085604316841996, 0.009139177302791153, 0.009191751407309061, 0.009243320884222412, 0.009293880097160322, 0.009343423520170318, 0.00939194573832249, 0.009439441448300992, 0.009485905458984029, 0.009531332692011101, 0.009575718182337876, 0.009619057078779164, 0.009661344644538864, 0.009702576257727958, 0.009742747411869234, 0.00978185371639035, 0.00981989089710335, 0.009856854796671895, 0.009892741375065664, 0.009927546710002122, 0.009961266997374775, 0.009993898551669387, 0.010025437806366347, 0.01005588131433098, 0.01008522574818994, 0.010113467900694973, 0.010140604685073598, 0.010166633135366204, 0.010191550406750467, 0.010215353775852252, 0.010238040641043068, 0.010259608522724598, 0.010280055063599757, 0.010299378028930195, 0.010317575306780569, 0.010334644908249447, 0.010350584967686661, 0.010365393742897078, 0.010379069615331337, 0.010391611090262378, 0.010403016796949039, 0.01041328548878585, 0.010422416043439223, 0.010430407462970154, 0.010437258873943282, 0.01044296952752239, 0.010447538799552137, 0.010450966190626462, 0.010453251326142981, 0.010454393956344087, 0.010454393956344087, 0.010453251326142981, 0.010450966190626462, 0.010447538799552137, 0.01044296952752239, 0.010437258873943282, 0.010430407462970154, 0.010422416043439223, 0.01041328548878585, 0.010403016796949039, 0.010391611090262378, 0.010379069615331337, 0.010365393742897078, 0.010350584967686661, 0.010334644908249447, 0.010317575306780569, 0.010299378028930195, 0.010280055063599757, 0.010259608522724598, 0.010238040641043068, 0.010215353775852252, 0.010191550406750467, 0.010166633135366204, 0.010140604685073598, 0.010113467900694973, 0.01008522574818994, 0.01005588131433098, 0.010025437806366347, 0.009993898551669387, 0.009961266997374775, 0.009927546710002122, 0.009892741375065664, 0.009856854796671895, 0.00981989089710335, 0.00978185371639035, 0.009742747411869234, 0.009702576257727958, 0.009661344644538864, 0.009619057078779164, 0.009575718182337876, 0.009531332692011101, 0.009485905458984029, 0.009439441448300992, 0.00939194573832249, 0.009343423520170318, 0.009293880097160322, 0.009243320884222412, 0.009191751407309061, 0.009139177302791153, 0.009085604316841996, 0.009031038304809223, 0.008975485230575371, 0.008918951165904783, 0.008861442289781127, 0.008802964887731923, 0.008743525351140993, 0.00868313017655067, 0.0086217859649506, 0.00855949942105765, 0.00849627735258214, 0.008432126669483946, 0.008367054383218206, 0.008301067605967079, 0.008234173549864183, 0.008166379526205846, 0.008097692944651201, 0.008028121312413226, 0.007957672233438363, 0.007886353407574561, 0.007814172629730136, 0.007741137789022482, 0.007667256867914979, 0.00759253794134379, 0.007516989175837821, 0.00744061882862333, 0.007363435246723478, 0.0072854468660452166, 0.007206662210456636, 0.007127089890856632, 0.00704673860423293, 0.0069656171327117265, 0.006883734342598478, 0.0068010991834066194, 0.006717720686882713, 0.006633607966017162, 0.006548770214047811, 0.006463216703456925, 0.00637695678495604, 0.0062899998864659286, 0.006202355512083302, 0.006114033241045046, 0.006025042726679271, 0.005935393695350927, 0.005845095945399296, 0.005754159346065033, 0.005662593836414576, 0.005570409424251355, 0.0054776161850220466, 0.0053842242607182395, 0.005290243858763116, 0.0051956852509013034, 0.0051005587720714925, 0.005004874819279474, 0.004908643850461695, 0.004811876383340857, 0.004714582994277626, 0.004616774317114607, 0.004518461042012614, 0.0044196539142855205, 0.004320363733222272, 0.004220601350909488, 0.004120377671042578, 0.004019703647735881, 0.003918590284326712, 0.0038170486321696183, 0.0037150897894316927, 0.0036127248998770794, 0.0035099651516518333, 0.003406821776058694, 0.003303306046330757, 0.0031994292764019384, 0.0030952028196669814, 0.0029906380677445906, 0.0028857464492311412, 0.002780539428452415, 0.002675028504212449, 0.0025692252085409055, 0.002463141105430975, 0.002356787789583519, 0.002250176885138956, 0.0021433200444213973, 0.0020362289466670158, 0.0019289152967680268, 0.0018213908240185923, 0.0017136672808561744, 0.0016057564416445678, 0.0014976701014677357, 0.0013894200749916825, 0.001281018195428937, 0.001172476313703897, 0.0010638062980465963, 0.0009550200344351969, 0.0008461294290789873, 0.0007371464159049574, 0.0006280829779724268, 0.0005189512134602094, 0.0004097635734033533, 0.00030053396541231163, 0.000191285544654933, 8.217779368318766e-05], } #legendre_points = { # 2: [-0.5773502691896257, 0.5773502691896257], # 3: [-0.7745966692414834, 0.0, 0.7745966692414834], # 4: [-0.8611363115940526, -0.33998104358485626, 0.33998104358485626, 0.8611363115940526], # 5: [-0.906179845938664, -0.5384693101056831, 0.0, 0.5384693101056831, 0.906179845938664], # 6: [-0.932469514203152, -0.6612093864662645, -0.23861918608319693, 0.23861918608319693, 0.6612093864662645, 0.932469514203152], # 7: [-0.9491079123427585, -0.7415311855993945, -0.4058451513773972, 0.0, 0.4058451513773972, 0.7415311855993945, 0.9491079123427585], # 8: [-0.9602898564975362, -0.7966664774136267, -0.525532409916329, -0.18343464249564978, 0.18343464249564978, 0.525532409916329, 0.7966664774136267, # 0.9602898564975362 # ], # 9: [-0.9681602395076261, -0.8360311073266358, -0.6133714327005904, -0.3242534234038089, 0.0, 0.3242534234038089, 0.6133714327005904, 0.8360311073266358, # 0.9681602395076261 # ], # 10: [-0.9739065285171717, -0.8650633666889845, -0.6794095682990244, -0.4333953941292472, -0.14887433898163122, 0.14887433898163122, 0.4333953941292472, # 0.6794095682990244, 0.8650633666889845, 0.9739065285171717 # ], # 11: [-0.978228658146057, -0.8870625997680953, -0.7301520055740494, -0.5190961292068118, -0.26954315595234496, 0.0, 0.26954315595234496, # 0.5190961292068118, 0.7301520055740494, 0.8870625997680953, 0.978228658146057 # ], # 12: [-0.9815606342467192, -0.9041172563704748, -0.7699026741943047, -0.5873179542866175, -0.3678314989981802, -0.1252334085114689, 0.1252334085114689, # 0.3678314989981802, 0.5873179542866175, 0.7699026741943047, 0.9041172563704748, 0.9815606342467192 # ], # 13: [-0.9841830547185881, -0.9175983992229779, -0.8015780907333099, -0.6423493394403402, -0.4484927510364468, -0.23045831595513477, 0.0, # 0.23045831595513477, 0.4484927510364468, 0.6423493394403402, 0.8015780907333099, 0.9175983992229779, 0.9841830547185881 # ], # 14: [-0.9862838086968123, -0.9284348836635735, -0.827201315069765, -0.6872929048116855, -0.5152486363581541, -0.31911236892788974, -0.10805494870734367, # 0.10805494870734367, 0.31911236892788974, 0.5152486363581541, 0.6872929048116855, 0.827201315069765, 0.9284348836635735, 0.9862838086968123 # ], # 15: [-0.9879925180204854, -0.937273392400706, -0.8482065834104272, -0.7244177313601701, -0.5709721726085388, -0.3941513470775634, -0.20119409399743451, # 0.0, 0.20119409399743451, 0.3941513470775634, 0.5709721726085388, 0.7244177313601701, 0.8482065834104272, 0.937273392400706, 0.9879925180204854 # ], # 16: [-0.9894009349916499, -0.9445750230732326, -0.8656312023878318, -0.755404408355003, -0.6178762444026438, -0.45801677765722737, -0.2816035507792589, - # 0.09501250983763745, 0.09501250983763745, 0.2816035507792589, 0.45801677765722737, 0.6178762444026438, 0.755404408355003, 0.8656312023878318, # 0.9445750230732326, 0.9894009349916499 # ], # 17: [-0.9905754753144174, -0.9506755217687678, -0.8802391537269859, -0.7815140038968014, -0.6576711592166908, -0.5126905370864769, -0.3512317634538763, - # 0.17848418149584785, 0.0, 0.17848418149584785, 0.3512317634538763, 0.5126905370864769, 0.6576711592166908, 0.7815140038968014, # 0.8802391537269859, 0.9506755217687678, 0.9905754753144174 # ], # 18: [-0.9915651684209309, -0.9558239495713978, -0.8926024664975557, -0.8037049589725231, -0.6916870430603532, -0.5597708310739475, -0.41175116146284263, - # 0.2518862256915055, -0.08477501304173529, 0.08477501304173529, 0.2518862256915055, 0.41175116146284263, 0.5597708310739475, 0.6916870430603532, # 0.8037049589725231, 0.8926024664975557, 0.9558239495713978, 0.9915651684209309 # ], # 19: [-0.9924068438435844, -0.96020815213483, -0.9031559036148179, -0.8227146565371428, -0.7209661773352294, -0.600545304661681, -0.46457074137596094, - # 0.31656409996362983, -0.1603586456402254, 0.0, 0.1603586456402254, 0.31656409996362983, 0.46457074137596094, 0.600545304661681, # 0.7209661773352294, 0.8227146565371428, 0.9031559036148179, 0.96020815213483, 0.9924068438435844 # ], ## 87: [-0.9996223504324843, -0.9980107215629823, -0.995113458137992, -0.9909336501109646, -0.9854766060684204, -0.9787493410038163, -0.9707605197472843, -0.9615204369360731, -0.9510410013893136, -0.9393357199742057, -0.9264196798912372, -0.9123095290941229, -0.8970234547675694, -0.8805811598520419, -0.8630038376297066, -0.8443141443970352, -0.8245361702557805, -0.8036954080581777, -0.7818187205454378, -0.7589343057213493, -0.7350716605053138, -0.7102615427114729, -0.6845359314028236, -0.6579279856713494, -0.6304720018972505, -0.6022033695423352, -0.5731585255345304, -0.5433749073022933, -0.5128909045194523, -0.4817458096226668, -0.44997976716528437, -0.4176337220728691, -0.3847493668670969, -0.3513690879260413, -0.31753591085011573, -0.28329344500409187, -0.24868582730667713, -0.2137576653401016, -0.17855397985304422, -0.14312014673100815, -0.10750183850894277, -0.07174496550149943, -0.035895616626800894, 0.0, 0.035895616626800894, 0.07174496550149943, 0.10750183850894277, 0.14312014673100815, 0.17855397985304422, 0.2137576653401016, 0.24868582730667713, 0.28329344500409187, 0.31753591085011573, 0.3513690879260413, 0.3847493668670969, 0.4176337220728691, 0.44997976716528437, 0.4817458096226668, 0.5128909045194523, 0.5433749073022933, 0.5731585255345304, 0.6022033695423352, 0.6304720018972505, 0.6579279856713494, 0.6845359314028236, 0.7102615427114729, 0.7350716605053138, 0.7589343057213493, 0.7818187205454378, 0.8036954080581777, 0.8245361702557805, 0.8443141443970352, 0.8630038376297066, 0.8805811598520419, 0.8970234547675694, 0.9123095290941229, 0.9264196798912372, 0.9393357199742057, 0.9510410013893136, 0.9615204369360731, 0.9707605197472843, 0.9787493410038163, 0.9854766060684204, 0.9909336501109646, 0.995113458137992, 0.9980107215629823, 0.9996223504324843], ## 300: [-0.9999679782184367, -0.9998312829844194, -0.9995853734488863, -0.9992302218632736, -0.998765860353095, -0.9981923380734435, -0.9975097171758793, -0.9967180720510002, -0.9958174891208242, -0.9948080667630034, -0.9936899152742682, -0.9924631568471006, -0.9911279255509899, -0.9896843673150097, -0.9881326399103604, -0.9864729129322606, -0.9847053677808957, -0.9828301976412724, -0.9808476074619004, -0.9787578139322587, -0.9765610454590226, -0.9742575421410374, -0.9718475557430343, -0.9693313496680812, -0.9667091989287728, -0.9639813901171553, -0.9611482213733913, -0.9582100023531646, -0.9551670541938287, -0.9520197094793031, -0.948768312203718, -0.9454132177338137, -0.9419547927700982, -0.9383934153067653, -0.93472947459038, -0.9309633710773336, -0.927095516390074, -0.9231263332721162, -0.9190562555418367, -0.9148857280450583, -0.9106152066064294, -0.9062451579796031, -0.9017760597962232, -0.8972084005137192, -0.8925426793619202, -0.88777940628849, -0.8829191019031915, -0.8779622974209857, -0.8729095346039718, -0.8677613657021737, -0.8625183533931814, -0.8571810707206516, -0.8517501010316768, -0.846226037913026, -0.8406094851262689, -0.8349010565417858, -0.8291013760716742, -0.8232110776015571, -0.817230804921301, -0.8111612116546523, -0.8050029611877982, -0.7987567265968609, -0.7924231905743324, -0.7860030453544584, -0.77949699263758, -0.7729057435134394, -0.7662300183834609, -0.7594705468820128, -0.752628067796661, -0.7457033289874226, -0.7386970873050261, -0.7316101085081913, -0.7244431671799336, -0.7171970466429053, -0.7098725388737805, -0.7024704444166952, -0.6949915722957499, -0.6874367399265869, -0.6798067730270485, -0.6721025055269295, -0.664324779476831, -0.6564744449561274, -0.6485523599800557, -0.6405593904059368, -0.6324964098385406, -0.6243642995346039, -0.6161639483065119, -0.6078962524251537, -0.5995621155219635, -0.5911624484901564, -0.5826981693851706, -0.5741702033243276, -0.5655794823857192, -0.5569269455063353, -0.5482135383794405, -0.5394402133512136, -0.530607929316659, -0.5217176516148028, -0.5127703519231844, -0.5037670081516558, -0.4947086043354983, -0.4855961305278722, -0.4764305826916062, -0.4672129625903424, -0.45794427767904683, -0.44862554099389806, -0.43925777104156566, -0.42984199168789095, -0.4203792320459818, -0.410870526363734, -0.40131691391079155, -0.3917194388649578, -0.38207915019807054, -0.37239710156135264, -0.36267435117025176, -0.35291196168878075, -0.34311100011337214, -0.333272537656259, -0.32339764962839507, -0.3134874153219266, -0.3035429178922298, -0.29356524423952524, -0.2835554848900835, -0.27351473387703407, -0.26344408862079155, -0.25334464980911087, -0.24321752127678611, -0.233063809885005, -0.22288462540037263, -0.21268108037361766, -0.2024542900179944, -0.19220537208739363, -0.18193544675417583, -0.17164563648674022, -0.16133706592684258, -0.15101086176667575, -0.140668152625726, -0.1303100689274186, -0.1199377427755664, -0.10955230783063467, -0.0991548991858358, -0.08874665324306741, -0.07832870758870747, -0.0679022008692799, -0.057468272667004265, -0.0470280633752434, -0.03658271407386212, -0.02613336640451104, -0.015681162445849026, -0.005227244588717748, 0.005227244588717748, 0.015681162445849026, 0.02613336640451104, 0.03658271407386212, 0.0470280633752434, 0.057468272667004265, 0.0679022008692799, 0.07832870758870747, 0.08874665324306741, 0.0991548991858358, 0.10955230783063467, 0.1199377427755664, 0.1303100689274186, 0.140668152625726, 0.15101086176667575, 0.16133706592684258, 0.17164563648674022, 0.18193544675417583, 0.19220537208739363, 0.2024542900179944, 0.21268108037361766, 0.22288462540037263, 0.233063809885005, 0.24321752127678611, 0.25334464980911087, 0.26344408862079155, 0.27351473387703407, 0.2835554848900835, 0.29356524423952524, 0.3035429178922298, 0.3134874153219266, 0.32339764962839507, 0.333272537656259, 0.34311100011337214, 0.35291196168878075, 0.36267435117025176, 0.37239710156135264, 0.38207915019807054, 0.3917194388649578, 0.40131691391079155, 0.410870526363734, 0.4203792320459818, 0.42984199168789095, 0.43925777104156566, 0.44862554099389806, 0.45794427767904683, 0.4672129625903424, 0.4764305826916062, 0.4855961305278722, 0.4947086043354983, 0.5037670081516558, 0.5127703519231844, 0.5217176516148028, 0.530607929316659, 0.5394402133512136, 0.5482135383794405, 0.5569269455063353, 0.5655794823857192, 0.5741702033243276, 0.5826981693851706, 0.5911624484901564, 0.5995621155219635, 0.6078962524251537, 0.6161639483065119, 0.6243642995346039, 0.6324964098385406, 0.6405593904059368, 0.6485523599800557, 0.6564744449561274, 0.664324779476831, 0.6721025055269295, 0.6798067730270485, 0.6874367399265869, 0.6949915722957499, 0.7024704444166952, 0.7098725388737805, 0.7171970466429053, 0.7244431671799336, 0.7316101085081913, 0.7386970873050261, 0.7457033289874226, 0.752628067796661, 0.7594705468820128, 0.7662300183834609, 0.7729057435134394, 0.77949699263758, 0.7860030453544584, 0.7924231905743324, 0.7987567265968609, 0.8050029611877982, 0.8111612116546523, 0.817230804921301, 0.8232110776015571, 0.8291013760716742, 0.8349010565417858, 0.8406094851262689, 0.846226037913026, 0.8517501010316768, 0.8571810707206516, 0.8625183533931814, 0.8677613657021737, 0.8729095346039718, 0.8779622974209857, 0.8829191019031915, 0.88777940628849, 0.8925426793619202, 0.8972084005137192, 0.9017760597962232, 0.9062451579796031, 0.9106152066064294, 0.9148857280450583, 0.9190562555418367, 0.9231263332721162, 0.927095516390074, 0.9309633710773336, 0.93472947459038, 0.9383934153067653, 0.9419547927700982, 0.9454132177338137, 0.948768312203718, 0.9520197094793031, 0.9551670541938287, 0.9582100023531646, 0.9611482213733913, 0.9639813901171553, 0.9667091989287728, 0.9693313496680812, 0.9718475557430343, 0.9742575421410374, 0.9765610454590226, 0.9787578139322587, 0.9808476074619004, 0.9828301976412724, 0.9847053677808957, 0.9864729129322606, 0.9881326399103604, 0.9896843673150097, 0.9911279255509899, 0.9924631568471006, 0.9936899152742682, 0.9948080667630034, 0.9958174891208242, 0.9967180720510002, 0.9975097171758793, 0.9981923380734435, 0.998765860353095, 0.9992302218632736, 0.9995853734488863, 0.9998312829844194, 0.9999679782184367], #} def jacobian(f, x0, scalar=True, perturbation=1e-9, zero_offset=1e-7, args=(), **kwargs): """def test_fun(x): # test case - 2 inputs, 3 outputs - should work fine x2 = x[0]*x[0] return np.array([x2*exp(x[1]), x2*sin(x[1]), x2*cos(x[1])]) def easy_fun(x): x = x[0] return 5*x*x - 3*x - 100 """ # For scalar - returns list, size of input variables # For vector - returns list of list - size of input variables * output variables # Could add backwards/complex, multiple evaluations, detection of poor condition # types and limits base = f(x0, *args, **kwargs) if not scalar: base = list(base) # copy the base point x = list(x0) nx = len(x0) gradient = [] for i in range(nx): delta = x0[i]*(perturbation) if delta == 0: delta = zero_offset x[i] += delta point = f(x, *args, **kwargs) if scalar: dy = (point - base)/delta gradient.append(dy) else: delta_inv = 1.0/delta dys = [delta_inv*(p - b) for p, b in zip(point, base)] gradient.append(dys) x[i] -= delta if not scalar: # Transpose to be in standard form return list(map(list, zip(*gradient))) return gradient def hessian(f, x0, scalar=True, perturbation=1e-9, zero_offset=1e-7, full=True, args=(), **kwargs): # Takes n**2/2 + 3*n/2 + 1 function evaluations! Can still be quite fast. # For scalar - returns list[list], size of input variables # For vector - returns list of list of list - size of input variables * input variables * output variables # Could add backwards/complex, multiple evaluations, detection of poor condition # types and limits, jacobian as output, fevals base = f(x0, *args, **kwargs) nx = len(x0) if not isinstance(base, (float, int, complex)): try: ny = len(base) except: ny = 1 else: ny = 1 deltas = [] for i in range(nx): delta = x0[i]*(perturbation) if delta == 0.0: delta = zero_offset deltas.append(delta) deltas_inv = [1.0/di for di in deltas] x_perturb = list(x0) fs_perturb_i = [] for i in range(nx): x_perturb[i] += deltas[i] f_perturb_i = f(x_perturb, *args, **kwargs) fs_perturb_i.append(f_perturb_i) x_perturb[i] -= deltas[i] if full: hessian = [[None]*nx for _ in range(nx)] else: hessian = [] for i in range(nx): if not full: row = [] f_perturb_i = fs_perturb_i[i] x_perturb[i] += deltas[i] for j in range(i+1): f_perturb_j = fs_perturb_i[j] x_perturb[j] += deltas[j] f_perturb_ij = f(x_perturb, *args, **kwargs) if scalar: dii0 = (f_perturb_i - base)*deltas_inv[i] dii1 = (f_perturb_ij - f_perturb_j)*deltas_inv[i] dij = (dii1 - dii0)*deltas_inv[j] else: # dii0s = [(fi - bi)*deltas_inv[i] for fi, bi in zip(f_perturb_i, base)] # dii1s = [(fij - fj)*deltas_inv[i] for fij, fj in zip(f_perturb_ij, f_perturb_j)] # dij = [(di1 - di0)*deltas_inv[j] for di1, di0 in zip(dii1s, dii0s)] # Saves a good amount of time dij = [((f_perturb_ij[m] - f_perturb_j[m]) - (f_perturb_i[m] - base[m]))*deltas_inv[j]*deltas_inv[i] for m in range(ny)] if not full: row.append(dij) else: hessian[i][j] = hessian[j][i] = dij x_perturb[j] -= deltas[j] if not full: hessian.append(row) x_perturb[i] -= deltas[i] return hessian def stable_poly_to_unstable(coeffs, low, high): if high != low: # Handle the case of no transformation, no limits my_poly = Polynomial([-0.5*(high + low)*2.0/(high - low), 2.0/(high - low)]) coeffs = horner(coeffs, my_poly).coef[::-1].tolist() return coeffs def horner_backwards(x, coeffs): return horner(coeffs, x) def exp_horner_backwards(x, coeffs): return exp(horner(coeffs, x)) def exp_horner_backwards_and_der(x, coeffs): poly_val, poly_der = horner_and_der(coeffs, x) val = exp(poly_val) der = poly_der*val return val, der def exp_horner_backwards_and_der2(x, coeffs): poly_val, poly_der, poly_der2 = horner_and_der2(coeffs, x) val = exp(poly_val) der = poly_der*val der2 = (poly_der*poly_der + poly_der2)*val return val, der, der2 def exp_horner_backwards_and_der3(x, coeffs): poly_val, poly_der, poly_der2, poly_der3 = horner_and_der3(coeffs, x) val = exp(poly_val) der = poly_der*val der2 = (poly_der*poly_der + poly_der2)*val der3 = (poly_der*poly_der*poly_der + 3.0*poly_der*poly_der2 + poly_der3)*val return val, der, der2, der3 def horner_backwards_ln_tau(T, Tc, coeffs): if T >= Tc: return 0.0 lntau = log(1.0 - T/Tc) return horner(coeffs, lntau) def horner_backwards_ln_tau_and_der(T, Tc, coeffs): if T >= Tc: return 0.0, 0.0 lntau = log(1.0 - T/Tc) val, poly_der = horner_and_der(coeffs, lntau) der = -poly_der/(Tc*(-T/Tc + 1)) return val, der def horner_backwards_ln_tau_and_der2(T, Tc, coeffs): if T >= Tc: return 0.0, 0.0, 0.0 lntau = log(1.0 - T/Tc) val, poly_der, poly_der2 = horner_and_der2(coeffs, lntau) der = -poly_der/(Tc*(-T/Tc + 1)) der2 = (-poly_der + poly_der2)/(Tc**2*(T/Tc - 1)**2) return val, der, der2 def horner_backwards_ln_tau_and_der3(T, Tc, coeffs): if T >= Tc: return 0.0, 0.0, 0.0, 00 lntau = log(1.0 - T/Tc) val, poly_der, poly_der2, poly_der3 = horner_and_der3(coeffs, lntau) der = -poly_der/(Tc*(-T/Tc + 1)) der2 = (-poly_der + poly_der2)/(Tc**2*(T/Tc - 1)**2) der3 = (2.0*poly_der - 3.0*poly_der2 + poly_der3)/(Tc**3*(T/Tc - 1)**3) return val, der, der2, der3 def exp_horner_backwards_ln_tau(T, Tc, coeffs): # This formulation has the nice property of being linear-linear when plotted # for surface tension if T >= Tc: return 0.0 # No matter what the polynomial term does to it, as tau goes to 1, x goes to a large negative value # So long as the polynomial has the right derivative at the end (and a reasonable constant) it will always converge to 0. lntau = log(1.0 - T/Tc) # Guarantee it is larger than 0 with the exp # This is a linear plot as well because both variables are transformed into a log basis. return exp(horner(coeffs, lntau)) def exp_horner_backwards_ln_tau_and_der(T, Tc, coeffs): if T >= Tc: return 0.0 tau = 1.0 - T/Tc lntau = log(tau) poly_val, poly_der_val = horner_and_der(coeffs, lntau) val = exp(poly_val) return val, -val*poly_der_val/(Tc*tau) def exp_horner_backwards_ln_tau_and_der2(T, Tc, coeffs): if T >= Tc: return 0.0 tau = 1.0 - T/Tc lntau = log(tau) poly_val, poly_val_der, poly_val_der2 = horner_and_der2(coeffs, lntau) val = exp(poly_val) der = -val*poly_val_der/(Tc*tau) der2 = (poly_val_der*poly_val_der - poly_val_der + poly_val_der2)*val/(Tc*Tc*(tau*tau)) return val, der, der2 def exp_poly_ln_tau_coeffs2(T, Tc, val, der): ''' from sympy import * T, Tc, T0, T1, T2, sigma0, sigma1, sigma2 = symbols('T, Tc, T0, T1, T2, sigma0, sigma1, sigma2') val, der = symbols('val, der') from sympy.abc import a, b, c from fluids.numerics import horner coeffs = [a, b] lntau = log(1 - T/Tc) sigma = exp(horner(coeffs, lntau)) d0 = diff(sigma, T) Eq0 = Eq(sigma,val) Eq1 = Eq(d0, der) s = solve([Eq0, Eq1], [a, b]) ''' c0 = der*(T - Tc)/val c1 = (-T*der*log((-T + Tc)/Tc) + Tc*der*log((-T + Tc)/Tc) + val*log(val))/val return [c0, c1] def exp_poly_ln_tau_coeffs3(T, Tc, val, der, der2): ''' from sympy import * T, Tc, T0, T1, T2, sigma0, sigma1, sigma2 = symbols('T, Tc, T0, T1, T2, sigma0, sigma1, sigma2') val, der, der2 = symbols('val, der, der2') from sympy.abc import a, b, c from fluids.numerics import horner coeffs = [a, b, c] lntau = log(1 - T/Tc) sigma = exp(horner(coeffs, lntau)) d0 = diff(sigma, T) Eq0 = Eq(sigma,val) Eq1 = Eq(d0, der) Eq2 = Eq(diff(d0, T), der2) # s = solve([Eq0, Eq1], [a, b]) s = solve([Eq0, Eq1, Eq2], [a, b, c]) ''' return [(-T**2*der**2 + T**2*der2*val + 2*T*Tc*der**2 - 2*T*Tc*der2*val + T*der*val - Tc**2*der**2 + Tc**2*der2*val - Tc*der*val)/(2*val**2), (T**2*der**2*log((-T + Tc)/Tc) - T**2*der2*val*log((-T + Tc)/Tc) - 2*T*Tc*der**2*log((-T + Tc)/Tc) + 2*T*Tc*der2*val*log((-T + Tc)/Tc) - T*der*val*log((-T + Tc)/Tc) + T*der*val + Tc**2*der**2*log((-T + Tc)/Tc) - Tc**2*der2*val*log((-T + Tc)/Tc) + Tc*der*val*log((-T + Tc)/Tc) - Tc*der*val)/val**2, (-T**2*der**2*log((-T + Tc)/Tc)**2 + T**2*der2*val*log((-T + Tc)/Tc)**2 + 2*T*Tc*der**2*log((-T + Tc)/Tc)**2 - 2*T*Tc*der2*val*log((-T + Tc)/Tc)**2 + T*der*val*log((-T + Tc)/Tc)**2 - 2*T*der*val*log((-T + Tc)/Tc) - Tc**2*der**2*log((-T + Tc)/Tc)**2 + Tc**2*der2*val*log((-T + Tc)/Tc)**2 - Tc*der*val*log((-T + Tc)/Tc)**2 + 2*Tc*der*val*log((-T + Tc)/Tc) + 2*val**2*log(val))/(2*val**2)] def horner_domain(x, coeffs, xmin, xmax): r'''Evaluates a polynomial defined by coefficienfs `coeffs` and domain (`xmin`, `xmax`) which maps the input variable into the window (-1, 1) where the polynomial can be evaluated most acccurately. The evaluation uses horner's method. Note that the coefficients are reversed compared to the common form; the first value is the coefficient of the highest-powered x term, and the last value in `coeffs` is the constant offset value. Parameters ---------- x : float Point at which to evaluate the polynomial, [-] coeffs : iterable[float] Coefficients of polynomial, [-] xmin : float Low value, [-] xmax : float High value, [-] Returns ------- val : float The evaluated value of the polynomial, [-] Notes ----- ''' range_inv = 1.0/(xmax - xmin) off = (-xmax - xmin)*range_inv scl = 2.0*range_inv x = off + scl*x tot = 0.0 for c in coeffs: tot = tot*x + c return tot def polynomial_offset_scale(xmin, xmax): range_inv = 1.0/(xmax - xmin) offset = (-xmax - xmin)*range_inv scale = 2.0*range_inv return offset, scale def horner_stable(x, coeffs, offset, scale): x = offset + scale*x tot = 0.0 for c in coeffs: tot = tot*x + c return tot def horner_stable_and_der(x, coeffs, offset, scale): x = offset + scale*x f = 0.0 der = 0.0 for a in coeffs: der = x*der + f f = x*f + a return (f, der*scale) def horner_stable_and_der2(x, coeffs, offset, scale): x = offset + scale*x f, der, der2 = 0.0, 0.0, 0.0 for a in coeffs: der2 = x*der2 + der der = x*der + f f = x*f + a return (f, der*scale, scale*scale*(der2 + der2)) def horner_stable_and_der3(x, coeffs, offset, scale): x = offset + scale*x f, der, der2, der3 = 0.0, 0.0, 0.0, 0.0 for a in coeffs: der3 = x*der3 + der2 der2 = x*der2 + der der = x*der + f f = x*f + a scale2 = scale*scale return (f, der*scale, scale2*(der2 + der2), scale2*scale*der3*6.0) def horner_stable_and_der4(x, coeffs, offset, scale): x = offset + scale*x f, der, der2, der3, der4 = 0.0, 0.0, 0.0, 0.0, 0.0 for a in coeffs: der4 = x*der4 + der3 der3 = x*der3 + der2 der2 = x*der2 + der der = x*der + f f = x*f + a scale2 = scale*scale return (f, der*scale, scale2*(der2 + der2), scale2*scale*der3*6.0, scale2*scale2*der4*24.0) def horner_stable_ln_tau(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0 lntau = log(1.0 - T/Tc) return horner_stable(lntau, coeffs, offset, scale) def horner_stable_ln_tau_and_der(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0, 0.0 lntau = log(1.0 - T/Tc) val, poly_der = horner_stable_and_der(lntau, coeffs, offset, scale) der = -poly_der/(Tc*(-T/Tc + 1)) return val, der def horner_stable_ln_tau_and_der2(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0, 0.0, 0.0 lntau = log(1.0 - T/Tc) val, poly_der, poly_der2 = horner_stable_and_der2(lntau, coeffs, offset, scale) der = -poly_der/(Tc*(-T/Tc + 1)) der2 = (-poly_der + poly_der2)/(Tc**2*(T/Tc - 1)**2) return val, der, der2 def horner_stable_ln_tau_and_der3(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0, 0.0, 0.0, 00 lntau = log(1.0 - T/Tc) val, poly_der, poly_der2, poly_der3 = horner_stable_and_der3(lntau, coeffs, offset, scale) der = -poly_der/(Tc*(-T/Tc + 1)) der2 = (-poly_der + poly_der2)/(Tc**2*(T/Tc - 1)**2) der3 = (2.0*poly_der - 3.0*poly_der2 + poly_der3)/(Tc**3*(T/Tc - 1)**3) return val, der, der2, der3 def exp_horner_stable(x, coeffs, offset, scale): return trunc_exp(horner_stable(x, coeffs, offset, scale)) def exp_horner_stable_and_der(x, coeffs, offset, scale): poly_val, poly_der = horner_stable_and_der(x, coeffs, offset, scale) val = exp(poly_val) der = poly_der*val return val, der def exp_horner_stable_and_der2(x, coeffs, offset, scale): poly_val, poly_der, poly_der2 = horner_stable_and_der2(x, coeffs, offset, scale) val = exp(poly_val) der = poly_der*val der2 = (poly_der*poly_der + poly_der2)*val return val, der, der2 def exp_horner_stable_and_der3(x, coeffs, offset, scale): poly_val, poly_der, poly_der2, poly_der3 = horner_stable_and_der3(x, coeffs, offset, scale) val = exp(poly_val) der = poly_der*val der2 = (poly_der*poly_der + poly_der2)*val der3 = (poly_der*poly_der*poly_der + 3.0*poly_der*poly_der2 + poly_der3)*val return val, der, der2, der3 def exp_horner_stable_ln_tau(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0 lntau = log(1.0 - T/Tc) return exp(horner_stable(lntau, coeffs, offset, scale)) def exp_horner_stable_ln_tau_and_der(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0, 0.0 tau = 1.0 - T/Tc lntau = log(tau) poly_val, poly_der_val = horner_stable_and_der(lntau, coeffs, offset, scale) val = exp(poly_val) return val, -val*poly_der_val/(Tc*tau) def exp_horner_stable_ln_tau_and_der2(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0, 0.0, 0.0 tau = 1.0 - T/Tc lntau = log(tau) poly_val, poly_val_der, poly_val_der2 = horner_stable_and_der2(lntau, coeffs, offset, scale) val = exp(poly_val) der = -val*poly_val_der/(Tc*tau) der2 = (poly_val_der*poly_val_der - poly_val_der + poly_val_der2)*val/(Tc*Tc*(tau*tau)) return val, der, der2 def exp_cheb(x, coeffs, offset, scale): y = chebval(x, coeffs, offset, scale) return trunc_exp(y) def exp_cheb_and_der(x, coeffs, der_coeffs, offset, scale): poly_val = chebval(x, coeffs, offset, scale) poly_der = chebval(x, der_coeffs, offset, scale) val = exp(poly_val) der = poly_der*val return val, der def exp_cheb_and_der2(x, coeffs, der_coeffs, der2_coeffs, offset, scale): poly_val = chebval(x, coeffs, offset, scale) poly_der = chebval(x, der_coeffs, offset, scale) poly_der2 = chebval(x, der2_coeffs, offset, scale) val = exp(poly_val) der = poly_der*val der2 = (poly_der*poly_der + poly_der2)*val return val, der, der2 def exp_cheb_and_der3(x, coeffs, der_coeffs, der2_coeffs, der3_coeffs, offset, scale): poly_val = chebval(x, coeffs, offset, scale) poly_der = chebval(x, der_coeffs, offset, scale) poly_der2 = chebval(x, der2_coeffs, offset, scale) poly_der3 = chebval(x, der3_coeffs, offset, scale) val = exp(poly_val) der = poly_der*val der2 = (poly_der*poly_der + poly_der2)*val der3 = (poly_der*poly_der*poly_der + 3.0*poly_der*poly_der2 + poly_der3)*val return val, der, der2, der3 def exp_cheb_ln_tau(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0 lntau = log(1.0 - T/Tc) y = chebval(lntau, coeffs, offset, scale) return trunc_exp(y) def exp_cheb_ln_tau_and_der(T, Tc, coeffs, der_coeffs, offset, scale): if T >= Tc: return 0.0, 0.0 tau = 1.0 - T/Tc lntau = log(tau) poly_val = chebval(lntau, coeffs, offset, scale) poly_der_val = chebval(lntau, der_coeffs, offset, scale) val = exp(poly_val) return val, -val*poly_der_val/(Tc*tau) def exp_cheb_ln_tau_and_der2(T, Tc, coeffs, der_coeffs, der2_coeffs, offset, scale): if T >= Tc: return 0.0, 0.0, 0.0 tau = 1.0 - T/Tc lntau = log(tau) poly_val = chebval(lntau, coeffs, offset, scale) poly_val_der = chebval(lntau, der_coeffs, offset, scale) poly_val_der2 = chebval(lntau, der2_coeffs, offset, scale) val = exp(poly_val) der = -val*poly_val_der/(Tc*tau) der2 = (poly_val_der*poly_val_der - poly_val_der + poly_val_der2)*val/(Tc*Tc*(tau*tau)) return val, der, der2 def chebval_ln_tau(T, Tc, coeffs, offset, scale): if T >= Tc: return 0.0 lntau = log(1.0 - T/Tc) poly_val = chebval(lntau, coeffs, offset, scale) return poly_val def chebval_ln_tau_and_der(T, Tc, coeffs, der_coeffs, offset, scale): if T >= Tc: return 0.0, 0.0 lntau = log(1.0 - T/Tc) val = chebval(lntau, coeffs, offset, scale) poly_der = chebval(lntau, der_coeffs, offset, scale) der = -poly_der/(Tc*(-T/Tc + 1)) return val, der def chebval_ln_tau_and_der2(T, Tc, coeffs, der_coeffs, der2_coeffs, offset, scale): if T >= Tc: return 0.0, 0.0, 0.0 lntau = log(1.0 - T/Tc) val = chebval(lntau, coeffs, offset, scale) poly_der = chebval(lntau, der_coeffs, offset, scale) poly_der2 = chebval(lntau, der2_coeffs, offset, scale) der = -poly_der/(Tc*(-T/Tc + 1)) der2 = (-poly_der + poly_der2)/(Tc**2*(T/Tc - 1)**2) return val, der, der2 def chebval_ln_tau_and_der3(T, Tc, coeffs, der_coeffs, der2_coeffs, der3_coeffs, offset, scale): if T >= Tc: return 0.0, 0.0, 0.0, 00 lntau = log(1.0 - T/Tc) val = chebval(lntau, coeffs, offset, scale) poly_der = chebval(lntau, der_coeffs, offset, scale) poly_der2 = chebval(lntau, der2_coeffs, offset, scale) poly_der3 = chebval(lntau, der3_coeffs, offset, scale) der = -poly_der/(Tc*(-T/Tc + 1)) der2 = (-poly_der + poly_der2)/(Tc**2*(T/Tc - 1)**2) der3 = (2.0*poly_der - 3.0*poly_der2 + poly_der3)/(Tc**3*(T/Tc - 1)**3) return val, der, der2, der3 def horner(coeffs, x): r'''Evaluates a polynomial defined by coefficienfs `coeffs` at a specified scalar `x` value, using the horner method. This is the most efficient formula to evaluate a polynomial (assuming non-zero coefficients for all terms). This has been added to the `fluids` library because of the need to frequently evaluate polynomials; and `NumPy`'s polyval is actually quite slow for scalar values. Note that the coefficients are reversed compared to the common form; the first value is the coefficient of the highest-powered x term, and the last value in `coeffs` is the constant offset value. Parameters ---------- coeffs : iterable[float] Coefficients of polynomial, [-] x : float Point at which to evaluate the polynomial, [-] Returns ------- val : float The evaluated value of the polynomial, [-] Notes ----- For maximum speed, provide a list of Python floats and `x` should also be of type `float` to avoid either `NumPy` types or slow python ints. Compare the speed with numpy via: >>> coeffs = np.random.uniform(0, 1, size=15) >>> coeffs_list = coeffs.tolist() %timeit np.polyval(coeffs, 10.0) `np.polyval` takes on the order of 15 us; `horner`, 1 us. Examples -------- >>> horner([1.0, 3.0], 2.0) 5.0 >>> horner([21.24288737657324, -31.326919865992743, 23.490607246508382, -14.318875366457021, 6.993092901276407, -2.6446094897570775, 0.7629439408284319, -0.16825320656035953, 0.02866101768198035, -0.0038190069303978003, 0.0004027586707189051, -3.394447111198843e-05, 2.302586717011523e-06, -1.2627393196517083e-07, 5.607585274731649e-09, -2.013760843818914e-10, 5.819957519561292e-12, -1.3414794055766234e-13, 2.430101267966631e-15, -3.381444175898971e-17, 3.4861255675373234e-19, -2.5070616549039004e-21, 1.122234904781319e-23, -2.3532795334141448e-26], 300.0) 1.9900667478569642e+58 References ---------- .. [1] "Horner’s Method." Wikipedia, October 6, 2018. https://en.wikipedia.org/w/index.php?title=Horner%27s_method&oldid=862709437. ''' tot = 0.0 for c in coeffs: tot = tot*x + c return tot def horner_and_der(coeffs, x): # Coefficients in same order as for horner f = 0.0 der = 0.0 for a in coeffs: der = x*der + f f = x*f + a return (f, der) def horner_and_der2(coeffs, x): # Coefficients in same order as for horner f, der, der2 = 0.0, 0.0, 0.0 for a in coeffs: der2 = x*der2 + der der = x*der + f f = x*f + a return (f, der, der2 + der2) def horner_and_der3(coeffs, x): # Coefficients in same order as for horner # Tested f, der, der2, der3 = 0.0, 0.0, 0.0, 0.0 for a in coeffs: der3 = x*der3 + der2 der2 = x*der2 + der der = x*der + f f = x*f + a return (f, der, der2 + der2, der3*6.0) def horner_and_der4(coeffs, x): # Coefficients in same order as for horner # Tested f, der, der2, der3, der4 = 0.0, 0.0, 0.0, 0.0, 0.0 for a in coeffs: der4 = x*der4 + der3 der3 = x*der3 + der2 der2 = x*der2 + der der = x*der + f f = x*f + a return (f, der, der2 + der2, der3*6.0, der4*24.0) def quadratic_from_points(x0, x1, x2, f0, f1, f2): ''' from sympy import * f, a, b, c, x, x0, x1, x2, f0, f1, f2 = symbols('f, a, b, c, x, x0, x1, x2, f0, f1, f2') func = a*x**2 + b*x + c Eq0 = Eq(func.subs(x, x0), f0) Eq1 = Eq(func.subs(x, x1), f1) Eq2 = Eq(func.subs(x, x2), f2) sln = solve([Eq0, Eq1, Eq2], [a, b, c]) cse([sln[a], sln[b], sln[c]], optimizations='basic', symbols=utilities.iterables.numbered_symbols(prefix='v')) ''' v0 = -x2 v1 = f0*(v0 + x1) v2 = f2*(x0 - x1) v3 = f1*(v0 + x0) v4 = x2*x2 v5 = x0*x0 v6 = x1*x1 v7 = 1.0/(v4*x0 + v5*x1 + v6*x2 - (v4*x1 + v5*x2 + v6*x0)) v8 = -v4 a = v7*(v1 + v2 - v3) b = -v7*(f0*(v6 + v8) - f1*(v5 + v8) + f2*(v5 - v6)) c = v7*(v1*x1*x2 + v2*x0*x1 - v3*x0*x2) return (a, b, c) def quadratic_from_f_ders(x, v, d1, d2): '''from sympy import * f, a, b, c, x, v, d1, d2 = symbols('f, a, b, c, x, v, d1, d2') f0 = a*x**2 + b*x + c f1 = diff(f0, x) f2 = diff(f0, x, 2) solve([Eq(f0, v), Eq(f1, d1), Eq(f2, d2)], [a, b, c]) ''' a = d2*0.5 b = d1 - d2*x c = -d1*x + d2*x*x*0.5 + v return [a, b, c] def is_poly_positive(poly, domain=None, rand_pts=10, j_tol=1e-12, root_perturb=1e-12): # Returns True if positive everywhere in the specified domain (or globally) if domain is None: # 1e-100 to 1e100 pts = logspace(-100, 100, rand_pts//2) pts += [-i for i in pts] else: pts = linspace(domain[0], domain[1], rand_pts) for p in pts: if horner(poly, p) < 0.0: return False roots = np.roots(poly) for root in roots: r = root.real if abs(root.imag/r) < j_tol: if (domain is not None) and (r < domain[0] or r > domain[1]): continue eps_high, eps_low = r*(1.0 + root_perturb), r*(1.0 - root_perturb) if horner(poly, eps_high) < 0: return False if horner(poly, eps_low) < 0: return False return True def is_poly_negative(poly, domain=None, rand_pts=10, j_tol=1e-12, root_perturb=1e-12): # Returns True if negative everywhere in the specified domain (or globally) poly = [-i for i in poly]# Changes the sign of all polynomial calculated values return is_poly_positive(poly, domain=domain, rand_pts=rand_pts, j_tol=j_tol, root_perturb=root_perturb) def min_max_ratios(actual, calculated): '''Given known and calculated data, compare the two and provide two numbers describing how far away from the known data the calculated data is. The numbers are the ratio of the lowest relative calc, and the highest relative calc. ''' min_ratio = max_ratio = 1.0 for i in range(len(actual)): if actual[i] == 0.0: r = 1.0 if calculated[i] == actual[i] else 10.0 else: r = calculated[i]/actual[i] if r < min_ratio: min_ratio = r elif r > max_ratio: max_ratio = r return min_ratio, max_ratio def std(data): '''Fast implementation of numpy's std function for 1d inputs only, with double precision only always. ''' tot = 0.0 N = len(data) for i in range(N): tot += data[i] mean = tot/N tot = 0.0 for i in range(N): v = (data[i] - mean) tot += v*v return sqrt(tot/N) def polyder(c, m=1, scl=1, axis=0): """not quite a copy of numpy's version because this was faster to implement.""" cnt = m if cnt == 0: return c n = len(c) if cnt >= n: c = c[:1]*0 else: for i in range(cnt): # normally only happens once n = n - 1 der = [0.0]*n for j in range(n, 0, -1): der[j - 1] = j*c[j] c = der return c def polyint(coeffs): """not quite a copy of numpy's version because this was faster to implement. Tried out a bunch of optimizations, and this hits a good balance between CPython and pypy speed.""" # return ([0.0] + [c/(i+1) for i, c in enumerate(coeffs[::-1])])[::-1] N = len(coeffs) out = [0.0]*(N+1) for i in range(N): out[i] = coeffs[i]/(N-i) return out def polyint_over_x(coeffs): N = len(coeffs) log_coef = coeffs[-1] Nm1 = N - 1 poly_terms = [0.0]*N for i in range(Nm1): poly_terms[i] = coeffs[i]/(Nm1-i) return poly_terms, log_coef # N = len(coeffs) # log_coef = coeffs[-1] # Nm1 = N - 1 # poly_terms = [coeffs[Nm1-i]/i for i in range(N-1, 0, -1)] # poly_terms.append(0.0) # return poly_terms, log_coef # coeffs = coeffs[::-1] # log_coef = coeffs[0] # poly_terms = [0.0] # for i in range(1, len(coeffs)): # poly_terms.append(coeffs[i]/i) # return list(reversed(poly_terms)), log_coef # def horner_log(coeffs, log_coeff, x): """Technically possible to save one addition of the last term of coeffs is removed but benchmarks said nothing was saved.""" tot = 0.0 for c in coeffs: tot = tot*x + c return tot + log_coeff*log(x) def fit_integral_linear_extrapolation(T1, T2, int_coeffs, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope): # Order T1, T2 so T2 is always larger for simplicity flip = T1 > T2 if flip: T1, T2 = T2, T1 tot = 0.0 if T1 < Tmin: T2_low = T2 if T2 < Tmin else Tmin x1 = Tmin_value - Tmin_slope*Tmin tot += T2_low*(0.5*Tmin_slope*T2_low + x1) - T1*(0.5*Tmin_slope*T1 + x1) if (Tmin <= T1 <= Tmax) or (Tmin <= T2 <= Tmax) or (T1 <= Tmin and T2 >= Tmax): T1_mid = T1 if T1 > Tmin else Tmin T2_mid = T2 if T2 < Tmax else Tmax tot += (horner(int_coeffs, T2_mid) - horner(int_coeffs, T1_mid)) if T2 > Tmax: T1_high = T1 if T1 > Tmax else Tmax x1 = Tmax_value - Tmax_slope*Tmax tot += T2*(0.5*Tmax_slope*T2 + x1) - T1_high*(0.5*Tmax_slope*T1_high + x1) if flip: return -tot return tot def poly_fit_integral_value(T, int_coeffs, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope): # Can still save 1 horner evaluation (all of them for height T), but will be VERY messy. if T < Tmin: x1 = Tmin_value - Tmin_slope*Tmin tot = T*(0.5*Tmin_slope*T + x1) return tot if (Tmin <= T <= Tmax): tot1 = horner(int_coeffs, T) - horner(int_coeffs, Tmin) x1 = Tmin_value - Tmin_slope*Tmin tot = Tmin*(0.5*Tmin_slope*Tmin + x1) return tot + tot1 else: x1 = Tmin_value - Tmin_slope*Tmin tot = Tmin*(0.5*Tmin_slope*Tmin + x1) tot1 = horner(int_coeffs, Tmax) - horner(int_coeffs, Tmin) x1 = Tmax_value - Tmax_slope*Tmax tot2 = T*(0.5*Tmax_slope*T + x1) - Tmax*(0.5*Tmax_slope*Tmax + x1) return tot1 + tot + tot2 def fit_integral_over_T_linear_extrapolation(T1, T2, T_int_T_coeffs, poly_fit_log_coeff, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope): # Order T1, T2 so T2 is always larger for simplicity flip = T1 > T2 if flip: T1, T2 = T2, T1 tot = 0.0 if T1 < Tmin: T2_low = T2 if T2 < Tmin else Tmin x1 = Tmin_value - Tmin_slope*Tmin tot += (Tmin_slope*T2_low + x1*log(T2_low)) - (Tmin_slope*T1 + x1*log(T1)) if (Tmin <= T1 <= Tmax) or (Tmin <= T2 <= Tmax) or (T1 <= Tmin and T2 >= Tmax): T1_mid = T1 if T1 > Tmin else Tmin T2_mid = T2 if T2 < Tmax else Tmax tot += (horner_log(T_int_T_coeffs, poly_fit_log_coeff, T2_mid) - horner_log(T_int_T_coeffs, poly_fit_log_coeff, T1_mid)) if T2 > Tmax: T1_high = T1 if T1 > Tmax else Tmax x1 = Tmax_value - Tmax_slope*Tmax tot += (Tmax_slope*T2 + x1*log(T2)) - (Tmax_slope*T1_high + x1*log(T1_high)) if flip: return -tot return tot def poly_fit_integral_over_T_value(T, T_int_T_coeffs, poly_fit_log_coeff, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope): if T < Tmin: x1 = Tmin_value - Tmin_slope*Tmin tot = (Tmin_slope*T + x1*log(T)) return tot if (Tmin <= T <= Tmax): tot1 = (horner_log(T_int_T_coeffs, poly_fit_log_coeff, T) - horner_log(T_int_T_coeffs, poly_fit_log_coeff, Tmin)) x1 = Tmin_value - Tmin_slope*Tmin tot = (Tmin_slope*Tmin + x1*log(Tmin)) return tot + tot1 else: x1 = Tmin_value - Tmin_slope*Tmin tot = (Tmin_slope*Tmin + x1*log(Tmin)) tot1 = (horner_log(T_int_T_coeffs, poly_fit_log_coeff, Tmax) - horner_log(T_int_T_coeffs, poly_fit_log_coeff, Tmin)) x2 = Tmax_value -Tmax*Tmax_slope tot2 = (-Tmax_slope*(Tmax - T) + x2*log(T) - x2*log(Tmax)) return tot1 + tot + tot2 def evaluate_linear_fits(data, x): calc = [] low_limits, high_limits, coeffs = data[0], data[3], data[6] for i in range(len(data[0])): if x < low_limits[i]: v = (x - low_limits[i])*data[1][i] + data[2][i] elif x > high_limits[i]: v = (x - high_limits[i])*data[4][i] + data[5][i] else: v = 0.0 for c in coeffs[i]: v = v*x + c # v = horner(coeffs[i], x) calc.append(v) return calc def evaluate_linear_fits_d(data, x): calc = [] low_limits, high_limits, dcoeffs = data[0], data[3], data[7] for i in range(len(data[0])): if x < low_limits[i]: dv = data[1][i] elif x > high_limits[i]: dv = data[4][i] else: dv = 0.0 for c in dcoeffs[i]: dv = dv*x + c calc.append(dv) return calc def evaluate_linear_fits_d2(data, x): calc = [] low_limits, high_limits, d2coeffs = data[0], data[3], data[8] for i in range(len(data[0])): d2v = 0.0 if low_limits[i] < x < high_limits[i]: for c in d2coeffs[i]: d2v = d2v*x + c calc.append(d2v) return calc def chebval(x, c, offset=0.0, scale=1.0): # Pure Python implementation of numpy.polynomial.chebyshev.chebval # This routine is faster in CPython as well as PyPy # Approxximately 2 adds and a multiply per coefficient x = offset + scale*x len_c = len(c) if len_c == 1: c0, c1 = c[0], 0.0 elif len_c == 2: c0, c1 = c[0], c[1] else: x2 = 2.0*x c0, c1 = c[-2], c[-1] for i in range(3, len_c + 1): c0_prev = c0 c0 = c[-i] - c1 c1 = c0_prev + c1*x2 return c0 + c1*x def chebder(c, m=1, scl=1.0): """not quite a copy of numpy's version because this was faster to implement. This does not evaluate the value of a cheb series at a point; it returns a new chebyshev seriese to be evaluated by chebval. """ c = list(c) cnt = int(m) if cnt == 0: return c n = len(c) if cnt >= n: c = [] else: for i in range(cnt): n = n - 1 if scl != 1.0: for j in range(len(c)): c[j] *= scl der = [0.0 for _ in range(n)] for j in range(n, 2, -1): der[j - 1] = (j + j)*c[j] c[j - 2] += (j*c[j])/(j - 2.0) if n > 1: der[1] = 4.0*c[2] der[0] = c[1] c = der return c def chebint(c, m=1, lbnd=0, scl=1): # k=[], is used by numpy to provide integration constants cnt = int(m) # the order of integration if cnt < 0: raise ValueError("Negative integration error not allowed") # if len(k) > cnt: # raise ValueError("Size of integration constants not consistent with order of integration") if cnt == 0: return c # k = list(k) + [0]*(cnt - len(k)) n = len(c) c2 = [0.0]*n # Make a copy of c for i in range(n): c2[i] = c[i] c = c2 for i in range(cnt): n = len(c) for o in range(n): c[o] *= scl if n == 1 and c[0] == 0: pass # c[0] += k[i] else: tmp = [0.0]*(n+1) tmp[1] = c[0] if n > 1: tmp[2] = c[1]/4 for j in range(2, n): cval = c[j] tmp[j + 1] = cval/(2.0*(j + 1.0)) tmp[j - 1] -= cval/(2.0*(j - 1.0)) # Scale is handled separately # tmp[0] += k[i] - chebval(lbnd, tmp) tmp[0] -= chebval(lbnd, tmp) c = tmp return c def binary_search(key, arr, size=None): if size is None: size = len(arr) imin = 0 imax = size if key > arr[size - 1]: return size while imin < imax: imid = imin + ((imax - imin) >> 1) if key >= arr[imid]: imin = imid + 1 else: imax = imid return imin - 1 def isclose(a, b, rel_tol=1e-9, abs_tol=0.0): """Pure python and therefore slow version of the standard library isclose. Works on older versions of python though! Hasn't been unit tested, but has been tested. manual unit testing: from math import isclose as isclose2 from random import uniform for i in range(10000000): a = uniform(-1, 1) b = uniform(-1, 1) rel_tol = uniform(0, 1) abs_tol = uniform(0, .001) ans1 = isclose(a, b, rel_tol, abs_tol) ans2 = isclose2(a, b, rel_tol=rel_tol, abs_tol=abs_tol) try: assert ans1 == ans2 except: print(a, b, rel_tol, abs_tol) """ if (rel_tol < 0.0 or abs_tol < 0.0 ): raise ValueError('Negative tolerances') if ((a.real == b.real) and (a.imag == b.imag)): return True if (isinf(a.real) or isinf(a.imag) or isinf(b.real) or isinf(b.imag)): return False diff = abs(a - b) return (((diff <= rel_tol*abs(b)) or (diff <= rel_tol*abs(a))) or (diff <= abs_tol)) try: from math import isclose except ImportError: pass def assert_close(a, b, rtol=1e-7, atol=0.0): if a is b: # Nice to handle None return True if __debug__: # Do not run these branches in -O, -OO mode try: try: assert isclose(a, b, rel_tol=rtol, abs_tol=atol) return except: assert cisclose(a, b, rel_tol=rtol, abs_tol=atol) return except: pass from numpy.testing import assert_allclose return assert_allclose(a, b, rtol=rtol, atol=atol) def assert_close1d(a, b, rtol=1e-7, atol=0.0): N = len(a) if N != len(b): raise ValueError("Variables are not the same length: %d, %d" %(N, len(b))) for i in range(N): assert_close(a[i], b[i], rtol=rtol, atol=atol) def assert_close2d(a, b, rtol=1e-7, atol=0.0): # N = len(a) # if N != len(b): # raise ValueError("Variables are not the same length: %d, %d" %(N, len(b))) # for i in range(N): # assert_close1d(a[i], b[i], rtol=rtol, atol=atol) N = len(a) if N != len(b): raise ValueError("Variables are not the same length: %d, %d" %(N, len(b))) if not __debug__: # Do not run these branches in -O, -OO mode from numpy.testing import assert_allclose return assert_allclose(a, b, rtol=rtol, atol=atol) for i in range(N): a0, b0 = a[i], b[i] N0 = len(a0) if N0 != len(b0): raise ValueError("Variables are not the same length: %d, %d" %(N0, len(b0))) for j in range(N0): # assert_close(a0[j], b0[j], rtol=rtol, atol=atol) good = True a1, b1 = a0[j], b0[j] if a1 is b1: # Nice to handle None pass else: try: try: good = isclose(a1, b1, rel_tol=rtol, abs_tol=atol) except: good = cisclose(a1, b1, rel_tol=rtol, abs_tol=atol) except: pass if not good: from numpy.testing import assert_allclose return assert_allclose(a1, b1, rtol=rtol, atol=atol) def assert_close3d(a, b, rtol=1e-7, atol=0.0): N = len(a) if N != len(b): raise ValueError("Variables are not the same length: %d, %d" %(N, len(b))) for i in range(N): assert_close2d(a[i], b[i], rtol=rtol, atol=atol) def assert_close4d(a, b, rtol=1e-7, atol=0.0): N = len(a) if N != len(b): raise ValueError("Variables are not the same length: %d, %d" %(N, len(b))) # for i in range(N): # assert_close3d(a[i], b[i], rtol=rtol, atol=atol) for i in range(N): a0, b0 = a[i], b[i] N0 = len(a0) if N0 != len(b0): raise ValueError("Variables are not the same length: %d, %d" %(N0, len(b0))) for j in range(N0): assert_close2d(a0[j], b0[j], rtol=rtol, atol=atol) def zeros(shape, dtype=float, fill_value=0.0): if dtype is float: zero = float(fill_value) elif dtype is int: zero = int(fill_value) else: raise ValueError("dtype is not supported") t = type(shape) if t is int: return [zero]*shape if t is tuple: l = len(shape) if l == 1: return [zero]*shape[0] elif l == 2: s1 = shape[1] return [[zero]*s1 for _ in range(shape[0])] elif l == 3: s1 = shape[1] s2 = shape[2] return [[[zero]*s2 for _ in range(s1)] for _ in range(shape[0])] elif l == 4: s1 = shape[1] s2 = shape[2] s3 = shape[3] return [[[[zero]*s3 for _ in range(s2)] for _ in range(s1)] for _ in range(shape[0])] raise ValueError("Shape not implemented") def full(shape, fill_value, dtype=float): return zeros(shape, dtype, fill_value) def interp(x, dx, dy, left=None, right=None, extrapolate=False): """One-dimensional linear interpolation routine inspired/ reimplemented from NumPy for extra speed for scalar values (and also numpy). Returns the one-dimensional piecewise linear interpolant to a function with a given value at discrete data-points. Parameters ---------- x : float X-coordinate of the interpolated values, [-] dx : list[float] X-coordinates of the data points, must be increasing, [-] dy : list[float] Y-coordinates of the data points; same length as `dx`, [-] left : float, optional Value to return for `x < dx[0]`, default is `dy[0]`, [-] right : float, optional Value to return for `x > dx[-1]`, default is `dy[-1]`, [-] extrapolate : bool, optional If True, for the cases of `left` and/or `right` not given, a linear extrapolation will be performed outside of bounds, [-] Returns ------- y : float The interpolated value, [-] Notes ----- This function is "unsafe" in that it assumes the x-coordinates are increasing. It also does not check for nan's, that `dx` and `dy` are the same length, and that `x` is scalar. Performance is 40-50% of that of NumPy under CPython. Examples -------- >>> interp(2.5, [1, 2, 3], [3, 2, 0]) 1.0 """ lendx = len(dx) j = binary_search(x, dx, lendx) if (j == -1): if left is not None: return left elif extrapolate: j = 0 return (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j])*(x - dx[j]) + dy[j] else: return dy[0] elif (j == lendx - 1): return dy[j] elif (j == lendx): if right is not None: return right elif extrapolate: j = -2 return (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j])*(x - dx[j]) + dy[j] else: return dy[-1] else: return (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j])*(x - dx[j]) + dy[j] def interp2d_linear(x, y, xs, ys, vals): # Same as RectBivariateSpline, s=0, kx=1, ky=1 (and better performance) if y < ys[0]: i0, i1 = 0, 1 y_dat = ys[i0], ys[i1] elif y > ys[-1]: i0, i1 = -2, -1 else: for i in range(len(ys)): if ys[i] >= y: i0, i1 = i-1, i break y_low, y_high = ys[i0], ys[i1] v_low = interp(x, xs, vals[i0], extrapolate=True) v_high = interp(x, xs, vals[i1], extrapolate=True) return v_low + (y-y_low)/(y_high-y_low)*(v_high-v_low) try: _array = np.array Polynomial = np.polynomial.Polynomial except: pass def implementation_optimize_tck(tck, force_numpy=False): """Converts 1-d or 2-d splines calculated with SciPy's `splrep` or. `bisplrep` to a format for fastest computation - lists in PyPy, and numpy arrays otherwise. Only implemented for 3 and 5 length `tck`s. """ if IS_PYPY_OR_SKIP_DEPENDENCIES and not force_numpy: return tuple(tck) else: size = len(tck) if size == 3: return (_array(tck[0]), _array(tck[1]), tck[2]) elif size == 5: return (_array(tck[0]), _array(tck[1]), _array(tck[2]), tck[3], tck[4]) else: raise NotImplementedError def tck_interp2d_linear(x, y, z, kx=1, ky=1): if kx != 1 or ky != 1: raise ValueError("Only linear formulations are currently implemented") # copy is not a method of lists in python 2 x = list(x) x.insert(0, x[0]) x.append(x[-1]) y = list(y) y.insert(0, y[0]) y.append(y[-1]) # c needs to be transposed, and made 1d c = [z[j][i] for i in range(len(z[0])) for j in range(len(z))] tck = [x, y, c, 1, 1] return implementation_optimize_tck(tck) def caching_decorator(f, full=False): from functools import wraps cache = {} info_cache = {} wraps = my_wraps() @wraps(f) def wrapper(x, *args, **kwargs): has_info = 'info' in kwargs if x in cache: if 'info' in kwargs: kwargs['info'][:] = info_cache[x] return cache[x] err = f(x, *args, **kwargs) cache[x] = err if has_info: info_cache[x] = list(kwargs['info']) return err if full: return wrapper, cache, info_cache return wrapper def translate_bound_func(func, bounds=None, low=None, high=None): if bounds is not None: low = [i[0] for i in bounds] high = [i[1] for i in bounds] def new_f(x, *args, **kwargs): """Function for a solver to call when using the bounded variables.""" x = [float(i) for i in x] for i in range(len(x)): x[i] = (low[i] + (high[i] - low[i])/(1.0 + exp(-x[i]))) # Return the actual results return func(x, *args, **kwargs) def translate_into(x): x = [float(i) for i in x] for i in range(len(x)): x[i] = -log((high[i] - x[i])/(x[i] - low[i])) return x def translate_outof(x): x = [float(i) for i in x] for i in range(len(x)): x[i] = (low[i] + (high[i] - low[i])/(1.0 + exp(-x[i]))) return x return new_f, translate_into, translate_outof def translate_bound_jac(jac, bounds=None, low=None, high=None): if bounds is not None: low = [i[0] for i in bounds] high = [i[1] for i in bounds] def new_j(x): x_base = [float(i) for i in x] N = len(x) for i in range(N): x_base[i] = (low[i] + (high[i] - low[i])/(1.0 + exp(-x[i]))) jac_base = jac(x_base) try: jac_base = [i for i in jac_base] for i in range(N): v = (high[i] - low[i])*exp(-x[i])*jac_base[i] v *= (1.0 + exp(-x[i]))**-2 jac_base[i] = v return jac_base except: raise NotImplementedError("Fail") def translate_into(x): x = [float(i) for i in x] for i in range(len(x)): x[i] = -log((high[i] - x[i])/(x[i] - low[i])) return x def translate_outof(x): x = [float(i) for i in x] for i in range(len(x)): x[i] = (low[i] + (high[i] - low[i])/(1.0 + exp(-x[i]))) return x return new_j, translate_into, translate_outof def translate_bound_f_jac(f, jac, bounds=None, low=None, high=None, inplace_jac=False, as_np=False): if bounds is not None: low = [i[0] for i in bounds] high = [i[1] for i in bounds] exp_terms = [0.0]*len(low) def new_f_j(x, *args): x_base = [i for i in x] N = len(x) for i in range(N): exp_terms[i] = ei = trunc_exp(-x[i]) x_base[i] = (low[i] + (high[i] - low[i])/(1.0 + ei)) if jac is True: f_base, jac_base = f(x_base, *args) else: f_base = f(x_base, *args) jac_base = jac(x_base, *args) try: if type(jac_base[0]) is list or (isinstance(jac_base, np.ndarray) and len(jac_base.shape) == 2): if not inplace_jac: jac_base = [[j for j in i] for i in jac_base] for i in range(len(jac_base)): for j in range(len(jac_base[i])): # Checked numerically t = (1.0 + exp_terms[j]) jac_base[i][j] = (high[j] - low[j])*exp_terms[j]*jac_base[i][j]/(t*t) else: if not inplace_jac: jac_base = [i for i in jac_base] for i in range(N): t = (1.0 + exp_terms[i]) jac_base[i] = (high[i] - low[i])*exp_terms[i]*jac_base[i]/(t*t) if as_np: jac_base = np.array(jac_base) return f_base, jac_base except: raise NotImplementedError("Fail") def translate_into(x): #x = [float(i) for i in x] for i in range(len(x)): x[i] = -trunc_log((high[i] - x[i])/(x[i] - low[i])) return x def translate_outof(x): #x = [float(i) for i in x] for i in range(len(x)): x[i] = (low[i] + (high[i] - low[i])/(1.0 + trunc_exp(-x[i]))) return x return new_f_j, translate_into, translate_outof class OscillationError(Exception): """Error raised when a derivative-based method is not converging.""" class UnconvergedError(Exception): """Error raised when maxiter has been reached in an optimization problem.""" def __repr__(self): return ('UnconvergedError("Failed to converge; maxiter (%d) reached, value=%g, error %g)"' %(self.maxiter, self.point, self.err)) def __init__(self, message, iterations=None, err=None, point=None): super(UnconvergedError, self).__init__(message) self.point = point self.iterations = iterations self.err = err class SamePointError(UnconvergedError): """Error raised when two trial points in a root finding problem have the same error.""" def __repr__(self): return 'TODO' def __init__(self, message, iterations=None, err=None, q1=None, p1=None, q0=None, p0=None): super(UnconvergedError, self).__init__(message) self.q1 = q1 self.p1 = p1 self.q0 = q0 self.p0 = p0 self.iterations = iterations self.err = err class NoSolutionError(Exception): """Error raised when detected that there is no actual solution to a problem.""" class NotBoundedError(Exception): """Error raised when a bisection type algorithm fails because its initial bounds do not bound the solution.""" class DiscontinuityError(Exception): """Error raised when a bisection type algorithm fails because there is a discontinuity.""" def damping_maintain_sign(x, step, damping=1.0, factor=0.5): """Damping function which will maintain the sign of the variable being manipulated. If the step puts it at the other sign, the distance between `x` and `step` will be shortened by the multiple of `factor`; i.e. if factor is `x`, the new value of `x` will be 0 exactly. The provided `damping` is applied as well. Parameters ---------- x : float Previous value in iteration, [-] step : float Change in `x`, [-] damping : float, optional The damping factor to be applied always, [-] factor : float, optional If the calculated step changes sign, this factor will be used instead of the step, [-] Returns ------- x_new : float The new value in the iteration, [-] Notes ----- Examples -------- >>> damping_maintain_sign(100, -200, factor=.5) 50.0 """ if isinstance(x, list): return [damping_maintain_sign(x[i], step[i], damping, factor) for i in range(len(x))] positive = x > 0.0 step_x = x + step if (positive and step_x < 0) or (not positive and step_x > 0.0): # print('damping') step = -factor*x return x + step*damping def make_damp_initial(steps=5, damping=1.0, *args): steps_holder = [steps] def damping_func(x, step, damping=damping, *args): if steps_holder[0] <= 0: # Do not dampen at all if isinstance(x, list): return [xi + dxi for xi, dxi in zip(x, step)] return x + step else: steps_holder[0] -= 1 if isinstance(x, list): return [xi + dxi*damping for xi, dxi in zip(x, step)] return x + step*damping return damping_func def make_max_step_initial(max_step, steps=5, *args): steps_holder = [steps] def damping_func(x, step, *args): if steps_holder[0] <= 0: # Do not dampen at all if isinstance(x, list): return [xi + dxi for xi, dxi in zip(x, step)] return x + step else: steps_holder[0] -= 1 if isinstance(x, list): next = [] for i in range(len(x)): the_step = step[i] if abs(the_step) > max_step: next.append(x[i] + copysign(max_step[i], the_step)) else: next.append(x[i] + the_step) return next the_step = step if abs(the_step) > max_step: return x + copysign(max_step, the_step) return x + the_step return damping_func def oscillation_checking_wrapper(f, minimum_progress=0.3, both_sides=False, full=True, good_err=None): checker = oscillation_checker(minimum_progress=minimum_progress, both_sides=both_sides, good_err=good_err) wraps = my_wraps() @wraps(f) def wrapper(x, *args, **kwargs): err_test = err = f(x, *args, **kwargs) if not isinstance(err, (float, int, complex)): err_test = err[0] try: oscillating = checker(x, err_test) except: oscillating = False # Zero division error probably if oscillating: raise OscillationError("Oscillating") return err if full: return wrapper, checker return wrapper class oscillation_checker(object): def __init__(self, minimum_progress=0.3, both_sides=False, good_err=None): self.minimum_progress = minimum_progress self.both_sides = both_sides self.xs_neg = [] self.xs_pos = [] self.ys_neg = [] self.ys_pos = [] # Provide a number that if the error is under this, no longer be able to return False # For example, near phase boundaries newton could be bisecting as it overshoots # each step, but is still converging fine self.good_err = good_err def clear(self): self.xs_neg = [] self.xs_pos = [] self.ys_neg = [] self.ys_pos = [] def is_solve_oscilating(self, x, y): if y == 0.0: return False xs_neg, xs_pos, ys_neg, ys_pos = self.xs_neg, self.xs_pos, self.ys_neg, self.ys_pos minimum_progress = self.minimum_progress if y < 0.0: xs_neg.append(x) ys_neg.append(y) else: xs_pos.append(x) ys_pos.append(y) if len(xs_pos) > 1 and len(xs_neg) > 1: if y < 0: dy_cur = y - ys_neg[-2] dy_other = ys_pos[-1] - ys_pos[-2] gain_neg = abs(dy_cur/y) gain_pos = abs(dy_other/ys_pos[-1]) else: dy_cur = y - ys_pos[-2] dy_other = ys_neg[-1] - ys_neg[-2] gain_pos = abs(dy_cur/y) gain_neg = abs(dy_other/ys_neg[-1]) # print(gain_pos, gain_neg, y) if self.both_sides: if gain_pos < minimum_progress and gain_neg < minimum_progress: if self.good_err is not None and min(abs(ys_neg[-1]), abs(ys_pos[-1])) < self.good_err: return False return True else: if gain_pos < minimum_progress or gain_neg < minimum_progress: if self.good_err is not None and min(abs(ys_neg[-1]), abs(ys_pos[-1])) < self.good_err: return False return True return False __call__ = is_solve_oscilating def best_bounding_bounds(low, high, f=None, xs_pos=None, ys_pos=None, xs_neg=None, ys_neg=None, fa=None, fb=None): r'''Given: 1) A presumed bracketing interval such as very far out limits on physical bounds 2) A history of a non-bounded search algorithm which did not converge Find the best bracketing points which get the algorithm as close to the solution as possible. Parameters ---------- low : float Low bracketing interval (`f` has opposite sign at `low` than `high`), [-] high : float High bracketing interval (`f` has opposite sign at `high` than `low`), [-] f : callable, optional 1D function to be solved, [-] xs_pos : list[float] Unsorted list of `x` values of points with positive `y` values previously evaluated, [-] ys_pos : list[float] Unsorted list of `y` values of points with positive `y` values previously evaluated, [-] xs_neg : list[float] Unsorted list of `x` values of points with negative `y` values previously evaluated, [-] ys_neg : list[float] Unsorted list of `y` values of points with negative `y` values previously evaluated, [-] fa : float, optional Value of function at `low`, used instead of recalculating if provided, [-] fb : float, optional Value of function at `high`, used instead of recalculating if provided, [-] Returns ------- low : float Low bracketing interval (`f` has opposite sign at `low` than `high`), [-] high : float High bracketing interval (`f` has opposite sign at `high` than `low`), [-] fa : float Value of function at `low`, [-] fb : float, optional Value of function at `high`, [-] Notes ----- Negative and/or positive history values can be omitted, but both the `x` and `y` lists should be skipped if so. More work could be done to handle better the case if the bounds not bracketing the root but the function history doing so. ''' if fa is None: fa = f(low) if fb is None: fb = f(high) if ys_pos: y_min_pos = min(ys_pos) x_min_pos = xs_pos[ys_pos.index(y_min_pos)] if fa > 0: if y_min_pos < fa: fa, low = y_min_pos, x_min_pos else: if y_min_pos < fb: fb, high = y_min_pos, x_min_pos if ys_neg: y_min_neg = max(ys_neg) x_min_neg = xs_neg[ys_neg.index(y_min_neg)] if fa < 0: if y_min_neg > fa: fa, low = y_min_neg, x_min_neg else: if y_min_pos > fb: fb, high = y_min_neg, x_min_neg if fa*fb > 0: raise ValueError("Bounds and previous history do not contain bracketing points") return low, high, fa, fb def bisect(f, a, b, args=(), xtol=1e-12, rtol=2.220446049250313e-16, maxiter=100, ytol=None): """Port of SciPy's C bisect routine.""" fa = f(a, *args) fb = f(b, *args) if fa*fb > 0.0: raise ValueError("f(a) and f(b) must have different signs") elif fa == 0.0: return a elif fb == 0.0: return b dm = b - a # iterations = 0.0 for i in range(maxiter): dm *= 0.5 xm = a + dm fm = f(xm, *args) if fm*fa >= 0.0: a = xm abs_dm = fabs(dm) if fm == 0.0: return xm elif ytol is not None: if (abs_dm < (xtol + rtol*abs_dm)) and abs(fm) < ytol: return xm elif (abs_dm < (xtol + rtol*abs_dm)): return xm # elif gap_detection: # dy_dx = abs((fm - fa)/(a-b)) # if dy_dx > dy_dx_limit: # raise DiscontinuityError("Discontinuity detected") raise UnconvergedError("Failed to converge after %d iterations" %maxiter) def ridder(f, a, b, args=(), xtol=_xtol, rtol=_rtol, maxiter=_iter, full_output=False, disp=True): a_abs, b_abs = fabs(a), fabs(b) tol = xtol + rtol*(a_abs if a_abs < b_abs else b_abs) fa = f(a, *args) fb = f(b, *args) if fa*fb > 0.0: raise ValueError("f(a) and f(b) must have different signs") elif fa == 0.0: return a elif fb == 0.0: return b for i in range(maxiter): dm = 0.5*(b - a) xm = a + dm fm = f(xm, *args) dn = copysign(1.0/sqrt(fm*fm - fa*fb), fb - fa)*fm*dm dn_abs, dm_abs_tol = fabs(dn), fabs(dm) - 0.5*tol xn = xm - copysign((dn_abs if dn_abs < dm_abs_tol else dm_abs_tol), dn) fn = f(xn, *args) if (fn*fm < 0.0): a = xn fa = fn b = xm fb = fm elif (fn*fa < 0.0): b = xn fb = fn else: a = xn fa = fn tol = xtol + rtol*xn if (fn == 0.0 or fabs(b - a) < tol): return xn raise UnconvergedError("Failed to converge after %d iterations" %maxiter) # numba: delete # raise UnconvergedError("Failed to converge") # numba: uncomment def brenth(f, xa, xb, args=(), xtol=1e-12, rtol=4.440892098500626e-16, maxiter=100, ytol=None, full_output=False, disp=True, q=False, fa=None, fb=None, kwargs={}): xpre = xa xcur = xb xblk = 0.0 fblk = 0.0 spre = 0.0 scur = 0.0 if fa is None: fpre = f(xpre, *args, **kwargs) else: fpre = fa if fb is None: fcur = f(xcur, *args, **kwargs) else: fcur = fb if fpre*fcur > 0.0: raise NotBoundedError("f(a) and f(b) must have different signs") elif fpre == 0.0: return xa elif fcur == 0.0: return xb for i in range(maxiter): if fpre*fcur < 0.0: xblk = xpre fblk = fpre spre = scur = xcur - xpre # Breaks a bunch of tests # if fpre == fcur: # raise UnconvergedError("Failed to converge - reached equal points after %d iterations" %i) if abs(fblk) < abs(fcur): xpre = xcur xcur = xblk xblk = xpre fpre = fcur fcur = fblk fblk = fpre delta = 0.5*(xtol + rtol*abs(xcur)) sbis = 0.5*(xblk - xcur) if ytol is not None: if fcur == 0.0 or (abs(sbis) < delta) and abs(fcur) < ytol: return xcur else: if fcur == 0.0 or (abs(sbis) < delta): return xcur if (abs(spre) > delta and abs(fcur) < abs(fpre)): if xpre == xblk: # interpolate stry = -fcur*(xcur - xpre)/(fcur - fpre) else: # extrapolate dpre = (fpre - fcur)/(xpre - xcur) dblk = (fblk - fcur)/(xblk - xcur) if q: stry = -fcur*(fblk*dblk - fpre*dpre)/(dblk*dpre*(fblk - fpre)) else: stry = -fcur*(fblk - fpre)/(fblk*dpre - fpre*dblk) if (abs(stry + stry) < min(abs(spre), 3.0*abs(sbis) - delta)): # accept step spre = scur scur = stry else: # bisect spre = sbis scur = sbis else: # bisect spre = sbis scur = sbis xpre = xcur fpre = fcur if abs(scur) > delta: xcur += scur else: if sbis > 0.0: xcur += delta else: xcur -= delta fcur = f(xcur, *args, **kwargs) raise UnconvergedError("Failed to converge after %d iterations" %maxiter) def secant(func, x0, args=(), maxiter=100, low=None, high=None, damping=1.0, xtol=1.48e-8, ytol=None, x1=None, require_eval=False, f0=None, f1=None, bisection=False, same_tol=1.0, kwargs={}, require_xtol=True): p0 = 1.0*x0 # Logic to take a small step to calculate the approximate derivative if x1 is not None: p1 = x1 else: if x0 >= 0.0: p1 = x0*1.0001 + 1e-4 else: p1 = x0*1.0001 - 1e-4 # May need to truncate p1 if low is not None and p1 < low: p1 = low if high is not None and p1 > high: p1 = high # Are we already converged on either point? Do not consider checking xtol # if so. if f0 is None: q0 = func(p0, *args, **kwargs) else: q0 = f0 if (ytol is not None and abs(q0) < ytol and not require_xtol) or q0 == 0.0: return p0 if f1 is None: q1 = func(p1, *args, **kwargs) else: q1 = f1 if (ytol is not None and abs(q1) < ytol and not require_xtol) or q1 == 0.0: return p1 if bisection: a, b = None, None if q1 < 0.0: a = p1 else: b = p1 if q0 < 0.0: a = p0 else: b = p0 for i in range(maxiter): # Calculate new point, and truncate if necessary if q1 != q0: p = p1 - q1*(p1 - p0)/(q1 - q0)*damping else: p = p1 if low is not None and p < low: p = low if high is not None and p > high: p = high # After computing new point if bisection and a is not None and b is not None: if not (a < p < b) or (b < p < a): p = 0.5*(a + b) # Check the exit conditions if ytol is not None and xtol is not None: # Meet both tolerance - new value is under ytol, and old value if abs(q1) < ytol and (not require_xtol or abs(p0 - p1) <= abs(xtol*p0)): # if abs(p0 - p1) <= abs(xtol*p0) and abs(q1) < ytol: if require_eval: return p1 return p elif xtol is not None: if abs(p0 - p1) <= abs(xtol*p0) and not (p0 == p1 and (p0 == low or p0 == high)): if require_eval: return p1 return p elif ytol is not None: if abs(q1) < ytol: if require_eval: return p1 return p # Check to quit after convergence check - may meet criteria if q1 == q0: # Are we close enough? Run the checks again if xtol is not None: xtol *= same_tol if ytol is not None: ytol *= same_tol if ytol is not None and xtol is not None: # Meet both tolerance - new value is under ytol, and old value if abs(p0 - p1) <= abs(xtol * p0) and abs(q1) < ytol: return p elif xtol is not None: if abs(p0 - p1) <= abs(xtol * p0) and not (p0 == p1 and (p0 == low or p0 == high)): return p elif ytol is not None: if abs(q1) < ytol: return p # Cannot proceed, raise an error raise SamePointError("Convergence failed - previous points are the same", q1=q1, p1=p1, q0=q0, p0=p0) # Swap the points around p0 = p1 q0 = q1 p1 = p q1 = func(p1, *args, **kwargs) if q1 == 0.0: return p1 if bisection: if q1 < 0.0: a = p1 else: b = p1 raise UnconvergedError("Failed to converge", iterations=i, point=p, err=q1) # start with 0.99 and try taking powers of two off of 1 line_search_factors = [] remove = 0.01 for i in range(7): line_search_factors.append(1-remove) remove*= 2 # once under 0.36, jump to 0.1 and go down by orders of magnitude to 1e-10 tmp = [] remove = 1e-10 for i in range(10): tmp.append(remove) remove *= 10 tmp.sort(reverse=True) line_search_factors.extend(tmp) line_search_factors_low_prec = line_search_factors[3:] newton_line_search_factors = [1.0] + line_search_factors newton_line_search_factors_disabled = [1.0] def one_sided_secant(f, x0, x_flat, args=tuple(), maxiter=100, xtol=1.48e-8, ytol=None, require_xtol=True, damping=1.0, x1=None, y_flat=None, max_quadratic_iter=7, low_prec_ls_iter=3, y0=None, y1=None, kwargs={}): if x1 is None: if x0 >= 0.0: x1 = x0*1.0001 + 1e-4 else: x1 = x0*1.0001 - 1e-4 if y0 is None: y0 = f(x0, *args, **kwargs) # print(f'Start point: x={x0}, y={y0}') if y1 is None: y1 = f(x1, *args, **kwargs) # print(f'Start point: x={x1}, y={y1}') flat_higher_than_x = x_flat > x0 # this is the flat value - all other y values that are flat must be this number # we only need to calculate it if not provided if y_flat is None: y_flat = f(x_flat, *args, **kwargs) if y0 == y_flat: raise ValueError("The initial y value is the same as the value in the zero-derivative region") if y0 == y1: raise ValueError("The initial points have the same value, cannot start the search") if y1 == y_flat: raise ValueError("The second point is also flat, cannot start the search") # store some other points for higher order steps x2 = x3 = y2 = y3 = None for i in range(maxiter): # # guess the new value based on the two points has_step = False if x2 is not None and i < max_quadratic_iter: try: # print('Points for quadratic', [x0, x1, x2], [y0, y1, y2]) a, b, c = quadratic_from_points(x0, x1, x2, y0, y1, y2) # print('Quadratic coefficients', (a, b, c)) sln = roots_quadratic(a, b, c) # print('Quadratic solutions', sln) if sln[0].imag == 0: if len(sln) == 2: a_step = (sln[0] - x1) another_step = (sln[1] - x1) # use the closest one if abs(a_step) < abs(another_step): # print('Using quadratic solution', sln[0]) step = a_step else: # print('Using quadratic solution', sln[1]) step = another_step else: # only got a single step # print('Using quadratic solution', sln[0]) step = (sln[0] - x1) has_step = True # print(f'Quadratic desired point: x={x1 + step}') except: # print('Quadratic had a numerical error') pass # secant can always work if not has_step: step = - y1*(x1 - x0)/(y1 - y0)*damping # print(f'Secant desired point: x={x1 + step}') force_line_search = x_flat <= (x1 + step) if flat_higher_than_x else x_flat >= (x1 + step) if not force_line_search: x = x1 + step y = f(x, *args, **kwargs) # print(f'Iteration: x={x}, y={y}, flat={y==y_flat}') if y == y_flat: x_flat = x else: pass # print(f'Secant desired point lower than flat point: x={x1 + step}, x_flat={x_flat}') if force_line_search or y == y_flat: # If the step is too big (overshoots the known flat point), recalculate it to be the limit if flat_higher_than_x: if x_flat <= (x0 + step): step = x_flat - x0 else: if x_flat >= (x0 + step): step = x_flat - x0 # Must use linesearch to find a x that gives a working y for line_search_factor in (line_search_factors if i > low_prec_ls_iter else line_search_factors_low_prec): x = x0 + step*line_search_factor y = f(x, *args, **kwargs) # print(f'Line search: x={x}, y={y}, flat={y==y_flat}, x_flat={x_flat}') if y != y_flat: break else: x_flat = x if y == y_flat: raise ValueError("The line search has finished and no point to continue with was found") # Check if we converged, but note this really requires a ytol as the xtol doesn't work well with lsearch if ytol is not None and xtol is not None: # Meet both tolerance - new value is under ytol, and old value if abs(y) < ytol and (abs(x - x0) <= abs(xtol*x)): return x elif xtol is not None: if abs(x - x0) <= abs(xtol*x): return x elif ytol is not None: if abs(y) < ytol: return x # change the points around, forgetting about one of the values x0, x1, x2, x3 = x, x0, x1, x2 y0, y1, y2, y3 = y, y0, y1, y2 raise UnconvergedError("Failed to converge", iterations=i, point=x, err=y) def halley_compat_numba(func, x, *args): a, b = func(x, *args) return a, b, 0.0 def newton(func, x0, fprime=None, args=(), tol=None, maxiter=100, fprime2=None, low=None, high=None, damping=1.0, ytol=None, xtol=1.48e-8, require_eval=False, damping_func=None, bisection=False, gap_detection=False, dy_dx_limit=1e100, max_bound_hits=4, kwargs={}): """Newton's method designed to be mostly compatible with SciPy's implementation, with a few features added and others now implemented. 1) No tracking of how many iterations have progressed. 2) No ability to return a RootResults object 3) No warnings on some cases of bad input (low tolerance, no iterations) 4) Ability to accept True for either fprime or fprime2, which means that they are included in the return value of func 5) No handling for inf or nan! 6) Special handling for functions which need to ensure an evaluation at the final point 7) Damping as a constant or a fraction 8) Ability to perform bisection, optionally specifying a maximum range 9) Ability to specify minimum and maximum iteration values 10) Ability to specify a tolerance in the `y` direction 11) Ability to pass in keyword arguments as well as positional arguments From scipy, with some modifications! https://github.com/scipy/scipy/blob/v1.1.0/scipy/optimize/zeros.py#L66-L206 """ if tol is not None: xtol = tol p0 = 1.0*x0 # p1 = p0 = 1.0*x0 # fval0 = None if bisection: a, b = None, None fa, fb = None, None fprime2_included = fprime2 == True fprime_included = fprime == True if fprime2_included: func2 = func hit_low, hit_high = 0, 0 for it in range(maxiter): # if fprime2_included: # numba: uncomment # fval, fder, fder2 = func(p0, *args) # numba: uncomment # else: # numba: uncomment # fval, fder = func(p0, *args) # numba: uncomment # fder2 = 0.0 # numba: uncomment if fprime2_included: # numba: DELETE fval, fder, fder2 = func2(p0, *args, **kwargs) # numba: DELETE elif fprime_included: # numba: DELETE fval, fder = func(p0, *args, **kwargs) elif fprime2 is not None: #numba: DELETE fval = func(p0, *args, **kwargs) #numba: DELETE fder = fprime(p0, *args, **kwargs) #numba: DELETE fder2 = fprime2(p0, *args, **kwargs) #numba: DELETE else: #numba: DELETE fval = func(p0, *args, **kwargs) #numba: DELETE fder = fprime(p0, *args, **kwargs) #numba: DELETE if fval == 0.0: return p0 # Cannot continue or already finished elif fder == 0.0: if ytol is None or abs(fval) < ytol: return p0 else: raise UnconvergedError("Derivative became zero; maxiter (%d) reached, value=%f " %(maxiter, p0)) if bisection: if fval < 0.0: a = p0 fa = fval else: b = p0 # b always has positive value of error fb = fval fder_inv = 1.0/fder # Compute the next point step = fval*fder_inv if damping_func is not None: if fprime2 is not None: step = step/(1.0 - 0.5*step*fder2*fder_inv) p = damping_func(p0, -step, damping) # variable, derivative, damping_factor elif fprime2 is None: p = p0 - step*damping else: p = p0 - step/(1.0 - 0.5*step*fder2*fder_inv)*damping if bisection and a is not None and b is not None: if (not (a < p < b) and not (b < p < a)): # if p < 0.0: # if p < a: # print('bisecting') p = 0.5*(a + b) # else: # if p > b: # p = 0.5*(a + b) if gap_detection: # Change in function value required to get goal in worst case dy_dx = abs((fa- fb)/(a-b)) if dy_dx > dy_dx_limit: #or dy_dx > abs(fder)*10: raise DiscontinuityError("Discontinuity detected") if low is not None and p < low: hit_low += 1 if p0 == low and hit_low > max_bound_hits: if abs(fval) < ytol: return low else: raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%f " % (maxiter, p)) # Stuck - not going to converge, hammering the boundary. Check ytol p = low if high is not None and p > high: hit_high += 1 if p0 == high and hit_high > max_bound_hits: if abs(fval) < ytol: return high else: raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%f " % (maxiter, p)) p = high # p0 is last point (fval at that point), p is new if ytol is not None and xtol is not None: # Meet both tolerance - new value is under ytol, and old value if abs(p - p0) < abs(xtol*p) and abs(fval) < ytol: if require_eval: return p0 return p elif xtol is not None: if abs(p - p0) < abs(xtol*p): if require_eval: return p0 return p elif ytol is not None: if abs(fval) < ytol: if require_eval: return p0 return p # fval0, fval1 = fval, fval0 # p0, p1 = p, p0 # need a history of fval also p0 = p # else: # return secant(func, x0, args=args, maxiter=maxiter, low=low, high=high, # damping=damping, # xtol=xtol, ytol=ytol, kwargs=kwargs) # raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%f " %(maxiter, p)) def halley(func, x0, args=(), maxiter=100, low=None, high=None, damping=1.0, ytol=None, xtol=1.48e-8, require_eval=False, damping_func=None, bisection=False, max_bound_hits=4, kwargs={}, max_2nd_ratio=1.5): p0 = 1.0*x0 if bisection: a, b = None, None fa, fb = None, None hit_low, hit_high = 0, 0 for it in range(maxiter): fval, fder, fder2 = func(p0, *args) if fval == 0.0: return p0 # Cannot continue or already finished elif fder == 0.0: if ytol is None or abs(fval) < ytol: return p0 else: raise UnconvergedError("Derivative became zero; maxiter (%d) reached, value=%f " %(maxiter, p0)) if bisection: if fval < 0.0: a = p0 fa = fval else: b = p0 # b always has positive value of error fb = fval fder_inv = 1.0/fder # Compute the next point step = fval*fder_inv skipped_halley = False if damping_func is not None: step = step/(1.0 - 0.5*step*fder2*fder_inv) p = damping_func(p0, -step, damping) else: step2 = step/(1.0 - 0.5*step*fder2*fder_inv)*damping if max_2nd_ratio is not None and abs(1.0-step2/step) > max_2nd_ratio: skipped_halley = True p = p0 - step else: p = p0 - step2 if bisection and a is not None and b is not None: if (not (a < p < b) and not (b < p < a)): p = 0.5*(a + b) if low is not None and p < low: hit_low += 1 if p0 == low and hit_low > max_bound_hits: if abs(fval) < ytol: return low else: raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%f " % (maxiter, p)) # Stuck - not going to converge, hammering the boundary. Check ytol p = low if high is not None and p > high: hit_high += 1 if p0 == high and hit_high > max_bound_hits: if abs(fval) < ytol: return high else: raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%f " % (maxiter, p)) p = high # p0 is last point (fval at that point), p is new if not skipped_halley: if ytol is not None and xtol is not None: # Meet both tolerance - new value is under ytol, and old value if abs(p - p0) < abs(xtol*p) and abs(fval) < ytol: if require_eval: return p0 return p elif xtol is not None: if abs(p - p0) < abs(xtol*p): if require_eval: return p0 return p elif ytol is not None: if abs(fval) < ytol: if require_eval: return p0 return p p0 = p raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%f " %(maxiter, p)) def newton_err(F): err = sum([abs(i) for i in F]) return err def basic_damping(x, dx, damping, *args): N = len(x) x2 = [0.0]*N for i in range(N): x2[i] = x[i] + dx[i]*damping return x2 def solve_2_direct(mat, vec): ab = mat[0] cd = mat[1] a, b = ab[0], ab[1] c, d = cd[0], cd[1] e, f = vec[0], vec[1] x0 = 1.0/(a*d - b*c) sln = [0.0]*2 sln[0] = x0*(-b*f + d*e) sln[1] = x0*(a*f - c*e) return sln def solve_3_direct(mat, vec): a, b, c = mat[0] d, e, f = mat[1] g, h, i = mat[2] j, k, l = vec x0 = a*e x1 = -b*d + x0 x2 = a*i - c*g x3 = a*f - c*d x4 = a*h - b*g x5 = x1*x2 - x3*x4 x6 = 1.0/x5 x7 = b*x3 - c*x1 x8 = 1.0/x1 x9 = d*x4 - g*x1 x10 = j*x8 x11 = a*k x12 = a*l ans = [0.0]*3 ans[0] = x6*(-k*x8*(b*x5 + x4*x7) + l*x7 + x10*(x0*x5 + x7*x9)/a) ans[1] = x6*(-x10*(d*x5 + x3*x9) + x11*x2 - x12*x3) ans[2] = x6*(j*x9 + x1*x12 - x11*x4) return ans def solve_4_direct(mat, vec): a, b, c, d = mat[0] e, f, g, h = mat[1] i, j, k, l = mat[2] m, n, o, p = mat[3] q, r, s, t = vec x0 = a*f x1 = -b*e + x0 x2 = x1*(a*k - c*i) x3 = a*g - c*e x4 = a*j - b*i x5 = x2 - x3*x4 x6 = a*h - d*e x7 = a*n - b*m x8 = x1*(a*p - d*m) - x6*x7 x9 = x1*(a*l - d*i) - x4*x6 x10 = x1*(a*o - c*m) - x3*x7 x11 = -x10*x9 + x5*x8 x12 = 1.0/x11 x13 = -b*x3 + c*x1 x14 = x13*x9 x15 = x5*(-b*x6 + d*x1) x16 = -x14 + x15 x17 = 1.0/x5 x18 = s*x17 x19 = x10*x4 - x5*x7 x20 = 1.0/x1 x21 = x17*x20 x22 = e*x4 - i*x1 x23 = x5*(e*x7 - m*x1) x24 = x10*x22 x25 = x23 - x24 x26 = q*x17 x27 = x20*x26 x28 = x3*x9 x29 = x5*x6 x30 = a*t x31 = -x28 + x29 x32 = a*r x33 = a*s*x1 x34 = x1*x30 x35 = -x23 + x24 ans = [0.0]*4 ans[0] = x12*(-r*x21*(-x11*(-b*x5 + x13*x4) + x16*x19) + t*(x14 - x15) - x18*(-x10*x16 + x11*x13) - x27*(-x11*(x0*x5 - x13*x22) + x16*x25)/a) ans[1] = x12*(-a*x18*(-x10*x31 + x11*x3) - x21*x32*(-x11*x2 + x19*x31) - x27*(x11*(e*x5 + x22*x3) + x25*x31) + x30*(x28 - x29)) ans[2] = x12*(-x17*x32*(x11*x4 + x19*x9) + x26*(x11*x22 + x35*x9) + x33*x8 - x34*x9) ans[3] = x12*(-q*x35 - x10*x33 + x19*x32 + x34*x5) return ans def newton_system(f, x0, jac, xtol=None, ytol=None, maxiter=100, damping=1.0, args=(), damping_func=None, solve_func=py_solve, line_search=False): # numba: delete # args=(), damping_func=None, solve_func=np.linalg.solve, line_search=False): # numba: uncomment jac_also = True if jac == True else False if jac_also: fcur, j = f(x0, *args) else: # numba: delete fcur = f(x0, *args) # numba: delete err0 = 0.0 for v in fcur: err0 += abs(v) if xtol is None and (ytol is not None and err0 < ytol): return x0, 0 else: x = x0 if not jac_also: # numba: delete j = jac(x, *args) # numba: delete factors = newton_line_search_factors if line_search else newton_line_search_factors_disabled iteration = 1 while iteration < maxiter: dx = solve_func(j, [-v for v in fcur]) # numba: delete # dx = solve_func(j, -fcur) # numba: uncomment for factor in factors: mult = factor*damping if damping_func is None: # xnew = x + dx*mult # numba: uncomment xnew = [xi + dxi*mult for xi, dxi in zip(x, dx)] # numba: delete else: xnew = damping_func(x, dx, damping*factor, *args) # try: if jac_also: fnew, jnew = f(xnew, *args) else: # numba: delete fnew = f(xnew, *args) # numba: delete # except: # print(f'Line search calculation with point failed') # continue err_new = 0.0 for v in fnew: err_new += abs(v) # print(f'Line search with error={err_new}, factor {mult}' ) if err_new < err0: break if line_search and err_new > err0: raise ValueError("Completed line search without reducing the objective function error, cannot proceed") fcur = fnew j = jnew err0 = err_new x = xnew iteration += 1 if xtol is not None: if (norm2(fcur) < xtol) and (ytol is None or err0 < ytol): break elif ytol is not None: if err0 < ytol: break if not jac_also: # numba: delete j = jac(x, *args) # numba: delete if xtol is not None and norm2(fcur) > xtol: raise UnconvergedError("Failed to converge") # raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%s" %(maxiter, x)) if ytol is not None: err0 = 0.0 for v in fcur: err0 += abs(v) if err0 > ytol: # raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%s" %(maxiter, x)) raise UnconvergedError("Failed to converge") return x, iteration def newton_minimize(f, x0, jac, hess, xtol=None, ytol=None, maxiter=100, damping=1.0, args=(), damping_func=None): jac_also = True if jac == True else False hess_also = True if hess == True else False def err(F): err = sum([abs(i) for i in F]) return err if hess_also: fcur, j, h = f(x0, *args) elif jac_also: fcur, j = f(x0, *args) h = hess(x0, *args) else: fcur = f(x0, *args) j = jac(x0, *args) h = hess(x0, *args) iter = 0 x = x0 while iter < maxiter: dx = py_solve(h, [-v for v in j]) if damping_func is None: x = [xi + dxi*damping for xi, dxi in zip(x, dx)] else: x = damping_func(x, dx, damping) if hess_also: fcur, j, h = f(x, *args) elif jac_also: fcur, j = f(x, *args) h = hess(x, *args) else: fcur = f(x, *args) j = jac(x, *args) h = hess(x, *args) iter += 1 if xtol is not None and norm2(j) < xtol: break if ytol is not None and err(j) < ytol: break if xtol is not None and norm2(j) > xtol: raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%s " %(maxiter, x)) if ytol is not None and err(j) > ytol: raise UnconvergedError("Failed to converge; maxiter (%d) reached, value=%s " %(maxiter, x)) return x, iter def broyden2(xs, fun, jac, xtol=1e-7, maxiter=100, jac_has_fun=False, skip_J=False, args=()): iter = 0 if skip_J: fcur = fun(xs, *args) N = len(fcur) J = eye(N) elif jac_has_fun: fcur, J = jac(xs, *args) J = inv(J) else: fcur = fun(xs, *args) J = inv(jac(xs, *args)) N = len(fcur) eqns = range(N) err = 0.0 for fi in fcur: err += abs(fi) while err > xtol and iter < maxiter: s = dot(J, fcur) xs = [xs[i] - s[i] for i in eqns] fnew = fun(xs, *args) z = [fnew[i] - fcur[i] for i in eqns] u = dot(J, z) d = [-i for i in s] dmu = [d[i]-u[i] for i in eqns] dmu_d = inner_product(dmu, d) den_inv = 1.0/inner_product(d, u) factor = den_inv*dmu_d J_delta = [[factor*j for j in row] for row in J] for i in eqns: for j in eqns: J[i][j] += J_delta[i][j] fcur = fnew iter += 1 err = 0.0 for fi in fcur: err += abs(fi) return xs, iter def normalize(values): r'''Simple function which normalizes a series of values to be from 0 to 1, and for their sum to add to 1. .. math:: x = \frac{x}{sum_i x_i} Parameters ---------- values : array-like array of values Returns ------- fractions : array-like Array of values from 0 to 1 Notes ----- Does not work on negative values, or handle the case where the sum is zero. Examples -------- >>> normalize([3, 2, 1]) [0.5, 0.3333333333333333, 0.16666666666666666] ''' tot_inv = 1.0/sum(values) return [i*tot_inv for i in values] # NOTE: the first value of each array is used; it is only the indexes that # are adjusted for fortran def func_35_splev(arg, t, l, l1, k2, nk1): # minus 1 index if arg >= t[l-1] or l1 == k2: arg, t, l, l1, nk1, leave = func_40_splev(arg, t, l, l1, nk1) # Always leaves here return arg, t, l, l1, k2, nk1 l1 = l l = l - 1 arg, t, l, l1, k2, nk1 = func_35_splev(arg, t, l, l1, k2, nk1) return arg, t, l, l1, k2, nk1 def func_40_splev(arg, t, l, l1, nk1): if arg < t[l1-1] or l == nk1: # minus 1 index return arg, t, l, l1, nk1, 1 l = l1 l1 = l + 1 arg, t, l, l1, nk1, leave = func_40_splev(arg, t, l, l1, nk1) return arg, t, l, l1, nk1, leave def py_splev(x, tck, ext=0, t=None, c=None, k=None): """Evaluate a B-spline using a pure-python port of FITPACK's splev. This is not fully featured in that it does not support calculating derivatives. Takes the knots and coefficients of a B-spline tuple, and returns the value of the smoothing polynomial. Parameters ---------- x : float or list[float] An point or array of points at which to calculate and return the value of the spline, [-] tck : 3-tuple Ssequence of length 3 returned by `splrep` containing the knots, coefficients, and degree of the spline, [-] ext : int, optional, default 0 If `x` is not within the range of the spline, this handles the calculation of the results. * For ext=0, extrapolate the value * For ext=1, return 0 as the value * For ext=2, raise a ValueError on that point and fail to return values * For ext=3, return the boundary value as the value Returns ------- y : list The array of calculated values; the spline function evaluated at the points in `x`, [-] Notes ----- The speed of this for a scalar value in CPython is approximately 15% slower than SciPy's FITPACK interface. In PyPy, this is 10-20 times faster than using it (benchmarked on PyPy 6). There could be more bugs in this port. """ e = ext if tck is not None: t, c, k = tck x = [x] # if isinstance(x, (float, int, complex)): # x = [x] m = 1#len(x) n = len(t) y = [] # output array k1 = k + 1 k2 = k1 + 1 nk1 = n - k1 tb = t[k1-1] # -1 to get right index te = t[nk1 ] # -1 to get right index; + 1 - 1 l = k1 l1 = l + 1 for i in range(0, m): # m is only 1 # i only used in loop for 1 arg = x[i] if arg < tb or arg > te: if e == 0: arg, t, l, l1, k2, nk1 = func_35_splev(arg, t, l, l1, k2, nk1) elif e == 1: y.append(0.0) continue elif e == 2: raise ValueError("X value not in domain; set `ext` to 0 to " "extrapolate") elif e == 3: if arg < tb: arg = tb else: arg = te arg, t, l, l1, k2, nk1 = func_35_splev(arg, t, l, l1, k2, nk1) else: arg, t, l, l1, k2, nk1 = func_35_splev(arg, t, l, l1, k2, nk1) # Local arrays used in fpbspl and to carry its result h = [0.0]*20 hh = [0.0]*19 fpbspl(t, n, k, arg, l, h, hh) sp = 0.0E0 ll = l - k1 for j in range(0, k1): ll = ll + 1 sp = sp + c[ll-1]*h[j] # -1 to get right index y.append(sp) return y[0] # if len(y) == 1: # return y[0] # return y def py_bisplev(x, y, tck, dx=0, dy=0): """Evaluate a bivariate B-spline or its derivatives. For scalars, returns a float; for other inputs, mimics the formats of SciPy's `bisplev`. Parameters ---------- x : float or list[float] x value (rank 1), [-] y : float or list[float] y value (rank 1), [-] tck : tuple(list, list, list, int, int) Tuple of knot locations, coefficients, and the degree of the spline, [tx, ty, c, kx, ky], [-] dx : int, optional Order of partial derivative with respect to `x`, [-] dy : int, optional Order of partial derivative with respect to `y`, [-] Returns ------- values : float or list[list[float]] Calculated values from spline or their derivatives; according to the same format as SciPy's `bisplev`, [-] Notes ----- Use `bisplrep` to generate the `tck` representation; there is no Python port of it. """ tx, ty, c, kx, ky = tck if isinstance(x, (float, int)): x = [x] if isinstance(y, (float, int)): y = [y] z = [[cy_bispev(tx, ty, c, kx, ky, [xi], [yi])[0] for yi in y] for xi in x] if len(x) == len(y) == 1: return z[0][0] return z def fpbspl(t, n, k, x, l, h, hh): """subroutine fpbspl evaluates the (k+1) non-zero b-splines of degree k at t(l) <= x < t(l+1) using the stable recurrence relation of de boor and cox. All arrays are 1d! Optimized the assignment and order and so on. """ h[0] = 1.0 for j in range(1, k + 1): hh[0:j] = h[0:j] h[0] = 0.0 for i in range(j): li = l+i f = hh[i]/(t[li] - t[li - j]) h[i] = h[i] + f*(t[li] - x) h[i + 1] = f*(x - t[li - j]) def init_w(t, k, x, lx, w): tb = t[k] n = len(t) m = len(x) h = [0]*6 hh = [0]*5 te = t[n - k - 1] l1 = k + 1 l2 = l1 + 1 for i in range(m): arg = x[i] if arg < tb: arg = tb if arg > te: arg = te while not (arg < t[l1] or l1 == (n - k - 1)): l1 = l2 l2 = l1 + 1 fpbspl(t, n, k, arg, l1, h, hh) lx[i] = l1 - k - 1 for j in range(k + 1): w[i][j] = h[j] def cy_bispev(tx, ty, c, kx, ky, x, y): """Possible optimization: Do not evaluate derivatives, ever.""" nx = len(tx) ny = len(ty) mx = len(x) my = len(y) kx1 = kx + 1 ky1 = ky + 1 nkx1 = nx - kx1 nky1 = ny - ky1 wx = [[0.0]*kx1]*mx wy = [[0.0]*ky1]*my lx = [0]*mx ly = [0]*my size_z = mx*my z = [0.0]*size_z init_w(tx, kx, x, lx, wx) init_w(ty, ky, y, ly, wy) for j in range(my): for i in range(mx): sp = 0.0 err = 0.0 for i1 in range(kx1): for j1 in range(ky1): l2 = lx[i]*nky1 + ly[j] + i1*nky1 + j1 a = c[l2]*wx[i][i1]*wy[j][j1] - err tmp = sp + a err = (tmp - sp) - a sp = tmp z[j*mx + i] += sp return z p_0_70 = array_if_needed([0.06184590404457956, -0.7460693871557973, 2.2435704485433376, -2.1944070385048526, 0.3382265629285811, 0.2791966558569478]) q_0_07 = array_if_needed([-0.005308735283483908, 0.1823421262956287, -1.2364596896290079, 2.9897802200092296, -2.9365321202088004, 1.0]) p_07_099 = array_if_needed([7543860.817140365, -10254250.429758755, -4186383.973408412, 7724476.972409749, -3130743.609030545, 600806.068543299, -62981.15051292659, 3696.7937385473397, -114.06795167646395, 1.4406337969700391]) q_07_099 = array_if_needed([-1262997.3422452002, 10684514.56076485, -16931658.916668657, 10275996.02842749, -3079141.9506451315, 511164.4690136096, -49254.56172495263, 2738.0399260270983, -81.36790509581284, 1.0]) p_099_1 = array_if_needed([8.548256176424551e+34, 1.8485781239087334e+35, -2.1706889553798647e+34, 8.318563643438321e+32, -1.559802348661511e+31, 1.698939241177209e+29, -1.180285031647229e+27, 5.531049937687143e+24, -1.8085903366375877e+22, 4.203276811951035e+19, -6.98211620300421e+16, 82281997048841.92, -67157299796.61345, 36084814.54808544, -11478.108105137717, 1.6370226052761176]) q_099_1 = array_if_needed([-1.9763570499484274e+35, 1.4813997374958851e+35, -1.4773854824041134e+34, 5.38853721252814e+32, -9.882387315028929e+30, 1.0635231532999732e+29, -7.334629044071992e+26, 3.420655574477631e+24, -1.1147787784365177e+22, 2.584530363912858e+19, -4.285376337404043e+16, 50430830490687.56, -41115254924.43107, 22072284.971253656, -7015.799744041691, 1.0]) def polylog2(x): r'''Simple function to calculate PolyLog(2, x) from ranges 0 <= x <= 1, with relative error guaranteed to be < 1E-7 from 0 to 0.99999. This is a Pade approximation, with three coefficient sets with splits at 0.7 and 0.99. An exception is raised if x is under 0 or above 1. Parameters ---------- x : float Value to evaluate PolyLog(2, x) T Returns ------- y : float Evaluated result Notes ----- Efficient (2-4 microseconds). No implementation of this function exists in SciPy. Derived with mpmath's pade approximation. Required for the entropy integral of :obj:`thermo.heat_capacity.Zabransky_quasi_polynomial`. Examples -------- >>> polylog2(0.5) 0.5822405264516294 ''' if 0 <= x <= 0.7: p = p_0_70 q = q_0_07 offset = 0.26 elif 0.7 < x <= 0.99: p = p_07_099 q = q_07_099 offset = 0.95 elif 0.99 < x <= 1: p = p_099_1 q = q_099_1 offset = 0.999 else: raise ValueError('Approximation is valid between 0 and 1 only.') x = x - offset return horner(p, x)/horner(q, x) def max_abs_error(data, calc): max_err = 0.0 N = len(data) for i in range(N): diff = abs(data[i] - calc[i]) if diff > max_err: max_err = diff return max_err def max_abs_rel_error(data, calc): max_err = 0.0 N = len(data) for i in range(N): diff = abs((data[i] - calc[i])/data[i]) if diff > max_err: max_err = diff return max_err def max_squared_error(data, calc): max_err = 0.0 N = len(data) for i in range(N): diff = abs(data[i] - calc[i]) diff *= diff if diff > max_err: max_err = diff return max_err def max_squared_rel_error(data, calc): max_err = 0.0 N = len(data) for i in range(N): diff = abs((data[i] - calc[i])/data[i]) diff *= diff if diff > max_err: max_err = diff return max_err def mean_abs_error(data, calc): mean_err = 0.0 N = len(data) for i in range(N): mean_err += abs(data[i] - calc[i]) return mean_err/N def mean_abs_rel_error(data, calc): mean_err = 0.0 N = len(data) for i in range(N): if data[i] == 0.0: if calc[i] == 0.0: pass else: mean_err += 1.0 else: mean_err += abs((data[i] - calc[i])/data[i]) return mean_err/N def mean_squared_error(data, calc): mean_err = 0.0 N = len(data) for i in range(N): err = abs(data[i] - calc[i]) mean_err += err*err return mean_err/N def mean_squared_rel_error(data, calc): mean_err = 0.0 N = len(data) for i in range(N): err = abs((data[i] - calc[i])/data[i]) mean_err += err*err return mean_err/N global sp_root sp_root = None def root(*args, **kwargs): global sp_root if sp_root is None: from scipy.optimize import root as sp_root return sp_root(*args, **kwargs) global sp_minimize sp_minimize = None def minimize(*args, **kwargs): global sp_minimize if sp_minimize is None: from scipy.optimize import minimize as sp_minimize return sp_minimize(*args, **kwargs) global sp_fsolve sp_fsolve = None def fsolve(*args, **kwargs): global sp_fsolve if sp_fsolve is None: from scipy.optimize import fsolve as sp_fsolve return sp_fsolve(*args, **kwargs) global sp_curve_fit sp_curve_fit = None def curve_fit(*args, **kwargs): global sp_curve_fit if sp_curve_fit is None: from scipy.optimize import curve_fit as sp_curve_fit return sp_curve_fit(*args, **kwargs) global sp_leastsq sp_leastsq = None def leastsq(*args, **kwargs): global sp_leastsq if sp_leastsq is None: from scipy.optimize import leastsq as sp_leastsq return sp_leastsq(*args, **kwargs) global sp_differential_evolution sp_differential_evolution = None def differential_evolution(*args, **kwargs): global sp_differential_evolution if sp_differential_evolution is None: from scipy.optimize import differential_evolution as sp_differential_evolution return sp_differential_evolution(*args, **kwargs) global sp_lmder sp_lmder = None def lmder(*args, **kwargs): global sp_lmder if sp_lmder is None: from scipy.optimize._minpack import _lmder as sp_lmder return sp_lmder(*args, **kwargs) global sp_lmfit sp_lmfit = None def lmfit(*args, **kwargs): global sp_lmfit if sp_lmfit is None: from scipy.optimize._minpack import _lmfit as sp_lmfit return sp_lmfit(*args, **kwargs) def fixed_quad_Gauss_Kronrod(f, a, b, k_points, k_weights, l_weights, args): ''' Note: This type of function cannot be fast in numba right now, the function call is too slow! https://numba.pydata.org/numba-doc/dev/user/faq.html#can-i-pass-a-function-as-an-argument-to-a-jitted-function ''' val_gauss_kronrod = val_gauss_legendre = 0.0 diff = b - a fact = 0.5*diff N = len(k_points) # center = N//2 # sp = [0.0]*N # fx_gk = [0.0]*N # fx_gl = [0.0]*center k = 0 for i in range(N): x0 = 0.5*(1.0 - k_points[i]) x1 = 0.5*(1.0 + k_points[i]) # sp[i] = x = a*x0 + b*x1 # fx_gk[i] = y = f(x, *args) val_gauss_kronrod += fact*y*k_weights[i] if i%2: # fx_gl[k] = y val_gauss_legendre += fact*y*l_weights[k] k += 1 return val_gauss_kronrod, val_gauss_kronrod-val_gauss_legendre def quad_adaptive(f, a, b, args=(), kronrod_points=array_if_needed(kronrod_points[10]), kronrod_weights=array_if_needed(kronrod_weights[10]), legendre_weights=array_if_needed(legendre_weights[10]), epsrel=1.49e-8, epsabs=1.49e-8, depth=0, points=None): # Disregard `points` for now area, err_abs = fixed_quad_Gauss_Kronrod(f, a, b, kronrod_points, kronrod_weights, legendre_weights, args) # Match behavior, documented at https://www.johndcook.com/blog/2012/03/20/scipy-integration/ #and with a good test case at https://www.johndcook.com/blog/2012/03/20/scipy-integration/ if (abs(err_abs) < epsabs or (area == 0.0 or abs(err_abs/area) < epsrel)) or depth > 6: # print((a, b), area, abs(err_abs), epsabs, abs(err_abs/area), epsrel, depth) return area, err_abs mid = a + (b-a)*0.5 area_A, err_abs_A = quad_adaptive(f, a, mid, args, kronrod_points, kronrod_weights, legendre_weights, epsrel=epsrel, epsabs=epsabs*0.5, depth=depth+1) area_B, err_abs_B = quad_adaptive(f, mid, b, args, kronrod_points, kronrod_weights, legendre_weights, epsrel=epsrel, epsabs=epsabs*0.5, depth=depth+1) return area_A + area_B, abs(err_abs_A) + abs(err_abs_B) global sp_quad sp_quad = None def lazy_quad(f, a, b, args=(), epsrel=1.49e-08, epsabs=1.49e-8, **kwargs): global sp_quad if not IS_PYPY: if sp_quad is None: from scipy.integrate import quad as sp_quad return sp_quad(f, a, b, args, epsrel=epsrel, epsabs=epsabs, **kwargs) else: return quad_adaptive(f, a, b, args=args, epsrel=epsrel, epsabs=epsabs) # n = 300 # return fixed_quad_Gauss_Kronrod(f, a, b, kronrod_points[n], kronrod_weights[n], legendre_weights[n], args) def _call_nquad(x, func, range_funcs, epsrel, epsabs, *args): return nquad(func, range_funcs, args=(x,) +args , epsrel=epsrel, epsabs=epsabs) def nquad(func, ranges, args=(), epsrel=1.48e-8, epsabs=1.48e-8): my_low, my_high = ranges[-1](*args) if len(ranges) == 1: return quad_adaptive(func, my_low, my_high, args=args, epsrel=epsrel, epsabs=epsabs) # return quad_adaptive(_call_nquad, my_low, my_high, args=(func, ranges[:-1], epsrel, epsabs) +args, epsrel=epsrel, epsabs=epsabs) def dblquad(func, a, b, hfun, gfun, args=(), epsrel=1.48e-12, epsabs=1.48e-15): """Nominally working, but trying to use it has exposed the expected bugs in `quad_adaptive`.""" def inner_func(y, *args): full_args = (y,)+args quad_fluids = quad_adaptive(func, hfun(y, *args), gfun(y, *args), args=full_args, epsrel=epsrel, epsabs=epsabs)[0] # from scipy.integrate import quad as quad2 # quad_sp = quad2(func, hfun(y), gfun(y), args=full_args, epsrel=epsrel, epsabs=epsabs)[0] # print(quad_fluids, quad_sp, hfun(y), gfun(y), full_args, ) return quad_fluids # return quad(func, hfun(y), gfun(y), args=(y,)+args, epsrel=epsrel, epsabs=epsabs)[0] return quad_adaptive(inner_func, a, b, args=args, epsrel=epsrel, epsabs=epsabs) if IS_PYPY: quad = quad_adaptive else: quad = lazy_quad fit_minimization_targets = {'MeanAbsErr': mean_abs_error, 'MeanRelErr': mean_abs_rel_error, 'MeanSquareErr': mean_squared_error, 'MeanSquareRelErr': mean_squared_rel_error, 'MaxAbsErr': max_abs_error, 'MaxRelErr': max_abs_rel_error, 'MaxSquareErr': max_squared_error, 'MaxSquareRelErr': max_squared_rel_error, } def py_factorial(n): if n < 0: raise ValueError("Positive values only") factorial = 1 for i in range(2, n + 1): factorial *= i return factorial try: from math import factorial except: factorial = py_factorial # interp, horner, derivative methods (and maybe newton?) should always be used. if not IS_PYPY: def splev(*args, **kwargs): from scipy.interpolate import splev return splev(*args, **kwargs) def bisplev(*args, **kwargs): from scipy.interpolate import bisplev return bisplev(*args, **kwargs) else: splev, bisplev = py_splev, py_bisplev # Try out mpmath for special functions anyway has_scipy = False if not SKIP_DEPENDENCIES and not is_micropython and not is_ironpython: has_scipy = True # try: # import scipy # has_scipy = True # except ImportError: # has_scipy = False else: has_scipy = False erf = None try: from math import erf except ImportError: # python 2.6 or other implementations? pass def _lambertw_err(x, y): return x*exp(x) - y def py_lambertw(y, k=0): """For x > 0, the is always only one real solution For -1/e < x < 0, two real solutions. Besides compatibility with scipy, the result should have a complex part because micropython doesn't support .real on floats """ # Works for real inputs only, two main branches if k == 0: # Branches dead at -1 # -1 is hard limit for real in this branch # 700 is safe upper limit for exp # Input should be between -1 and +BIGNUMBER return brenth(_lambertw_err, -1.0, 700.0, (y,)) + 0.0j elif k == -1: # Input should be between 0 and -1/e # not a big input range! return brenth(_lambertw_err, -700.0, -1.0, (y,)) + 0.0j else: raise ValueError("Other branches not supported") #has_scipy = False if IS_PYPY: lambertw = py_lambertw if has_scipy: if not IS_PYPY: def lambertw(*args, **kwargs): from scipy.special import lambertw return lambertw(*args, **kwargs) def ellipe(m): from scipy.special import ellipe return ellipe(m) def gammaincc(*args, **kwargs): from scipy.special import gammaincc return gammaincc(*args, **kwargs) def i1(*args, **kwargs): from scipy.special import i1 return i1(*args, **kwargs) def i0(*args, **kwargs): from scipy.special import i0 return i0(*args, **kwargs) def k1(*args, **kwargs): from scipy.special import k1 return k1(*args, **kwargs) def k0(*args, **kwargs): from scipy.special import k0 return k0(*args, **kwargs) def iv(*args, **kwargs): from scipy.special import iv return iv(*args, **kwargs) def hyp2f1(*args, **kwargs): from scipy.special import hyp2f1 return hyp2f1(*args, **kwargs) def ellipkinc(phi, m): from scipy.special import ellipkinc return ellipkinc(phi, m) def ellipeinc(phi, m): from scipy.special import ellipeinc return ellipeinc(phi, m) if erf is None: def erf(*args, **kwargs): from scipy.special import erf return erf(*args, **kwargs) # from scipy.special import lambertw, ellipe, gammaincc # fluids # from scipy.special import i1, i0, k1, k0, iv # ht # from scipy.special import hyp2f1 # if erf is None: # from scipy.special import erf else: lambertw = py_lambertw # scipy is not available... fall back to mpmath as a Pure-Python implementation # However, lazy load all functions # from mpmath import lambertw # Same branches as scipy, supports .real # Figured out this definition from test_precompute_gammainc.py in scipy if erf is None: def erf(*args, **kwargs): import mpmath return mpmath.erf(*args, **kwargs) def gammaincc(a, x): import mpmath return mpmath.gammainc(a, a=x, regularized=True) def ellipe(*args, **kwargs): import mpmath return mpmath.ellipe(*args, **kwargs) def gammaincc(a, x): import mpmath return mpmath.gammainc(a, a=x, regularized=True) def iv(*args, **kwargs): import mpmath return mpmath.besseli(*args, **kwargs) def hyp2f1(*args, **kwargs): import mpmath return mpmath.hyp2f1(*args, **kwargs) def iv(*args, **kwargs): import mpmath return mpmath.besseli(*args, **kwargs) def i1(x): import mpmath return mpmath.mpmath.besseli(1, x) def i0(x): import mpmath return mpmath.mpmath.besseli(0, x) def k1(x): import mpmath return mpmath.mpmath.besselk(1, x) def k0(x): import mpmath return mpmath.mpmath.besselk(0, x) def ellipkinc(phi, m): import mpmath return mpmath.mpmath.ellipf(phi, m) def ellipeinc(phi, m): import mpmath return mpmath.ellipe.ellipeinc(phi, m) try: if FORCE_PYPY: lambertw = py_lambertw from scipy.special import ellipe, gammaincc, gamma, ellipkinc, ellipeinc # fluids from scipy.special import i1, i0, k1, k0, iv # ht from scipy.special import hyp2f1 if erf is None: from scipy.special import erf except: pass