# -*- 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, isnan, 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_product, norm2, matrix_vector_dot, eye, array_as_tridiagonals, tridiagonals_as_array, transpose, solve_tridiagonal, subset_matrix, argsort1d, shape, sort_paired_lists, stack_vectors, matrix_multiply, transpose, eye, inv, sum_matrix_rows, gelsd) from fluids.numerics.special import (py_hypot, py_cacos, py_catan, py_catanh, trunc_exp, trunc_log, cbrt, factorial, comb) from fluids.numerics.polynomial_roots import (roots_quadratic, roots_quartic, roots_cubic_a1, roots_cubic_a2, roots_cubic) from fluids.numerics.polynomial_evaluation import (horner, horner_and_der, horner_and_der2, horner_and_der3, horner_and_der4, horner_backwards, exp_horner_backwards, exp_horner_backwards_and_der, exp_horner_backwards_and_der2, exp_horner_backwards_and_der3, horner_backwards_ln_tau, horner_backwards_ln_tau_and_der, horner_backwards_ln_tau_and_der2, horner_backwards_ln_tau_and_der3, exp_horner_backwards_ln_tau, exp_horner_backwards_ln_tau_and_der, exp_horner_backwards_ln_tau_and_der2, horner_domain, horner_stable, horner_stable_and_der, horner_stable_and_der2, horner_stable_and_der3, horner_stable_and_der4, horner_stable_ln_tau, horner_stable_ln_tau_and_der, horner_stable_ln_tau_and_der2, horner_stable_ln_tau_and_der3, exp_horner_stable, exp_horner_stable_and_der, exp_horner_stable_and_der2, exp_horner_stable_and_der3, exp_horner_stable_ln_tau, exp_horner_stable_ln_tau_and_der, exp_horner_stable_ln_tau_and_der2, horner_log, horner_stable_log, ) from fluids.numerics.polynomial_utils import (polyint, polyint_over_x, polyder, quadratic_from_points, quadratic_from_f_ders, deflate_cubic_real_roots, exp_poly_ln_tau_coeffs3, exp_poly_ln_tau_coeffs2, polynomial_offset_scale, stable_poly_to_unstable, polyint_stable, polyint_over_x_stable, poly_convert) __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', 'exp_poly_ln_tau_coeffs3', 'exp_poly_ln_tau_coeffs2', 'implementation_optimize_tck', 'tck_interp2d_linear', 'bisect', 'ridder', 'brenth', 'newton', 'secant', 'halley', 'one_sided_secant', 'trunc_exp_numpy', 'trunc_log_numpy', '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', 'comb', 'SolverInterface', 'multivariate_solvers', 'jacobian_methods', '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', 'homotopy_solver', 'horner_stable_log', 'is_increasing', 'fixed_point_anderson', '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', 'fixed_point_to_residual', 'residual_to_fixed_point', 'sort_paired_lists', # 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', 'is_monotonic', 'sort_nelder_mead_points_numba', 'sort_nelder_mead_points_python', 'bounds_clip_naive', 'nelder_mead', 'cbrt', 'polyint_stable', 'polyint_over_x_stable', 'argsort1d', 'poly_convert', ] 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 try: from math import cbrt except: def cbrt(x): return x**(1.0/3.0) _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 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 central_diff_weights(points, divisions=1): # Check the cache key = (divisions, points) try: return central_diff_weights_precomputed[key] except: pass 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[key] = 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. """ try: weights = central_diff_weights_precomputed[(n, order)] except: if order < n + 1: raise ValueError if order % 2 == 0: raise ValueError weights = central_diff_weights(order, n) ho = order >> 1 if n == 1: denominator = 1.0/dx else: denominator = 1.0 for i in range(n): denominator *= dx denominator = 1.0/denominator if scalar: if upper_limit is not None: max_x = x0 + (order - 1 - ho)*dx if max_x > upper_limit: x0 -= (max_x - x0) if lower_limit is not None: min_x = x0 + -ho*dx if min_x < lower_limit: x0 += (x0 - min_x) # failed optimization strategy - plus changes some answers # const = x0 - dx*ho # tot += weights[k]*func(const + dx*k, *args, **kwargs) 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=(), base=None, **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 if base is None: 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, base=None, 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 if base is None: 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): tmp = x_perturb[i] x_perturb[i] += deltas[i] f_perturb_i = f(x_perturb, *args, **kwargs) fs_perturb_i.append(f_perturb_i) x_perturb[i] = tmp # specifically set the same value back - floating point prec can't add and subtract to get the exact same value always 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] tmp_i = x_perturb[i] x_perturb[i] += deltas[i] for j in range(i+1): f_perturb_j = fs_perturb_i[j] tmp_j = x_perturb[j] if i == j: x_perturb[j] = tmp_i + 2.0*deltas[j] else: x_perturb[j] += deltas[j] f_perturb_ij = f(x_perturb, *args, **kwargs) if scalar: # dii0 = (f_perturb_i - base) # dii1 = (f_perturb_ij - f_perturb_j) # dij = (dii1 - dii0)*deltas_inv[j]*deltas_inv[i] # This exact equation helps avoid errors when the equations are independent dij = ((f_perturb_ij+base) - (f_perturb_i+f_perturb_j) )*deltas_inv[j]*deltas_inv[i] 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] = tmp_j if not full: hessian.append(row) x_perturb[i] = tmp_i return hessian 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 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 is_monotonic(points): """Checks if the input sequence of points is monotonic. A sequence is monotonic if it is entirely increasing or ecreasing. This function iterates through the provided sequence to determine if it meets this criterion. Parameters ---------- points : list or array_like A list of numerical values representing a sequence of points, [-] Returns ------- monotonic : bool `True` if the sequence is monotonic, `False` if the sequence contains both increases and decreases. Examples -------- >>> is_monotonic([1, 2, 2, 3]) True >>> is_monotonic([3, 2, 1]) True >>> is_monotonic([1, 3, 2, 4]) False >>> is_monotonic([1e5, 1e7, 1e10]) True Notes ----- The function considers a single-element sequence as monotonic. It also treats sequences with equal adjacent values as monotonic (e.g., [1, 2, 2, 3]). """ N = len(points) if N == 1: return True a_point_increased = False a_point_decreased = False for i in range(N-1): if points[i+1] > points[i]: a_point_increased = True if points[i+1] < points[i]: a_point_decreased = True return not (a_point_increased and a_point_decreased) def is_increasing(points): """Check if a given sequence of points is entirely increasing. Points must be larger, not equal, by this criteria Parameters ---------- points : list or array_like A list of numerical values representing a sequence of points, [-] Returns ------- increasing : bool `True` if the sequence is entirely increasing (each element is strictly greater than its predecessor), `False` otherwise, [-] Examples -------- >>> is_increasing([1, 2, 3, 4]) True >>> is_increasing([1, 1, 2, 3]) False >>> is_increasing([4, 3, 2, 1]) False Notes ----- Empty sequences or length 1 sequences are considered increasing. """ N = len(points) if N < 2: return True previous = points[0] for i in range(1, N): current = points[i] if current <= previous: return False previous = current return True 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 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 trunc_exp_numpy(x, trunc=1.7976931348622732e+308, out=None): if out is not None: out.fill(trunc) else: out = np.full_like(x, trunc) mask = np.where(x< 709.782712893384) out[mask] = np.exp(x[mask]) return out def trunc_log_numpy(x, trunc=-744.4400719213812, out=None): if out is not None: out.fill(trunc) else: out = np.full_like(x, trunc) if np.any(x < 0.0): raise ValueError("math domain error") mask = np.where(x > 0.0) out[mask] = np.log(x[mask]) return out 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: return c[0] if len_c == 2: return c[0] + c[1] * x 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. """ if m == 0: return c c = list(c) n = len(c) if m >= n: return [] der = [0.0] * (n - 1) # Pre-allocate der with maximum possible size for _ in range(m): n = n - 1 if scl != 1.0: for j in range(n+1): c[j] *= scl 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 = der, c # Swap c and der for _ in range(len(c) - n): c.pop() return c def chebint(c, m=1, lbnd=0, scl=1): cnt = int(m) # the order of integration if cnt < 0: raise ValueError("Negative integration error not allowed") if cnt == 0: return c n0 = len(c) if n0 == 1 and c[0] == 0.0: return c nmax = n0 + cnt # Maximum size after integrations c0 = [0.0] * nmax # Pre-allocate c0 and tmp with maximum needed size tmp = [0.0] * nmax # Copy initial coefficients into c0 c0[:n0] = c[:] n = n0 # Current length of c0 for _ in range(cnt): # Scale c0 if necessary if scl != 1: for i in range(n): c0[i] *= scl if n == 1 and c0[0] == 0: # Special case where c0 is zero pass else: tmp[1] = c0[0] if n > 1: tmp[2] = 0.25 * c0[1] for j in range(2, n): cval = c0[j] tmp[j + 1] = cval / (2.0 * (j + 1)) tmp[j - 1] -= cval / (2.0 * (j - 1)) val = chebval(lbnd, tmp) tmp[0] -= val # Swap c0 and tmp for next iteration c0, tmp = tmp, c0 n += 1 # Degree increases by 1 after integration return c0 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)) # Voodoo to handle the fact numba doesn't have isclose try: IS_NUMBA except: 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 + trunc_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] = -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, 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 + trunc_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])*trunc_exp(-x[i])*jac_base[i] v *= (1.0 + trunc_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] = -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_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 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, 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 or isnan(fpre) or isnan(fcur): 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) factors_growing_positive = [3.0**i for i in range(1, 20)] factors_shrinking_positive = [3.0**(-i) for i in range(1, 20)] factors_growing_negative = [-v for v in factors_growing_positive] factors_shrinking_negative = [-v for v in factors_shrinking_positive] factors_nearby_increasing = [v for v in linspace(1, 3, 20)] factors_nearby_decreasing = [1.0/v for v in factors_nearby_increasing] secant_bisection_factors = [] for i in range(len(factors_growing_positive)): for l in (factors_growing_positive, factors_shrinking_positive, factors_growing_negative, factors_shrinking_negative, factors_nearby_increasing, factors_nearby_decreasing): if l[i] != 1.0: secant_bisection_factors.append(l[i]) 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, additional_guesses=False, additional_guess_scale=1.0): 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 isnan(q0) or isinf(q0): raise UnconvergedError("Bad function evaluation on first starting point", iterations=0, point=p0, err=q0) 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 isnan(q1) or isinf(q1): raise UnconvergedError("Bad function evaluation on second starting point", iterations=0, point=p1, err=q1) if q1 == q0 and additional_guesses: # we cannot proceed when the guessed point has the same value; need to increase the search space # try to guess again with a larger change, up to 16 orders of magnitude larger for guess_factor in secant_bisection_factors: # print('Slope was zero with initial guess, expanding search...') p1 = x0*guess_factor if x0 != 0.0 else guess_factor if low is not None and p1 < low: p1 = low if high is not None and p1 > high: p1 = high q1 = func(p1, *args, **kwargs) # For the check, we will require a significant difference in result magnitude if not (isnan(q1) or isinf(q1)) and abs(q1/q0-1.0) > 1e-10: break did_additional_guesses = False 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): # print('Bisecting') 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 if bisection and a is not None and b is not None: p = 0.5*(a + b) else: # Search increasingly away from the initial guess for something with a different sign if additional_guesses and not did_additional_guesses: # print('Converged to points with a slope of zero, attempting to bisect the problem...') for guess_factor in secant_bisection_factors: p = x0*guess_factor*additional_guess_scale if (low is not None and p < low) or (high is not None and p > high): continue temp_q = func(p, *args, **kwargs) # print(f'Did not bisect the problem at x={p}') if not (isnan(q1) or isinf(q1)) and temp_q*q0 < 0.0: # print('Bisected the problem, restarting with new point') break did_additional_guesses = True else: # 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 (isnan(q1) or isinf(q1)): raise UnconvergedError("Bad function evaluation", iterations=i, point=p1, err=q1) if q1 == 0.0: return p1 if bisection: if q1 < 0.0: a = p1 else: b = p1 # print(f"f({p1}) = {q1}, high={a}, low={b}") 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 sln[0] != sln[1]: a_step = (sln[0].real - x1) another_step = (sln[1].real - 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].real - 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=(), maxiter=100, fprime2=None, low=None, high=None, damping=1.0, ytol=None, xtol=1.48e-8, require_eval=False, require_xtol=True, damping_func=None, bisection=False, additional_guesses=False, 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 """ p0 = 1.0*x0 p1 = None p2 = None p3 = None # Yes triple loops occur fval0 = None fval1 = None fval2 = None fval3 = 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 it = 0 did_additional_guesses = False while it < 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, fcurv = 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 fcurv = fprime2(p0, *args, **kwargs) #numba: DELETE else: #numba: DELETE fval = func(p0, *args, **kwargs) #numba: DELETE fder = fprime(p0, *args, **kwargs) #numba: DELETE bad_fval = isnan(fval) or isinf(fval) bad_fder = fder == 0.0 or isnan(fder) or isinf(fder) # print(f'Guess {p0}, value {fval}, derivative {fder} on iteration {it}') stuck = False if p1 is not None and p2 is not None: # for numba stuck = p0 == p2 and p1 == p0# or (fval == fval0 and fval == fval1) if not stuck and p3 is not None: stuck = p3 == p0# or fval == fval3# triple loop if fval == 0.0: # print('Completed search, fval is zer0') return p0 # Cannot continue or already finished elif bad_fder or bad_fval or stuck: if not bad_fval and (ytol is not None and abs(fval) < ytol): return p0 if stuck and (not additional_guesses or did_additional_guesses): raise UnconvergedError("Failed to converge; next guess is same as current guess on iteration %d, guess=%g, value=%g"%(it, p0, fval)) if it == 0 and not additional_guesses: if bad_fval: raise UnconvergedError("Cannot continue - math error in function value on iteration %d, guess=%f, value=%f"%(it, p0, fval)) else: raise UnconvergedError("Derivative became zero on iteration %d, guess=%f " %(it, p0)) require_bracket = False if additional_guesses and ((it == 0 or stuck) or (fval == fval0)): # fval == fval0 checks if we are flat and a line search is pointless points = secant_bisection_factors reinitializing = True did_additional_guesses = True require_bracket = fval == fval0 and bisection# The derivative led us astray, try the other way else: points = line_search_factors reinitializing = False # print(f'Guess {p0} could not progress value={fval} on iteration {it} , expanding search...') recovered = False for guess_factor in points: it += 1 if reinitializing: guess = x0*guess_factor if x0 != 0.0 else guess_factor else: guess = p1 + (p0 - p1)*guess_factor if low is not None and guess < low: guess = low if high is not None and guess > high: guess = high if guess == p0: continue if fprime2_included: # numba: DELETE fval, fder, fcurv = func2(guess, *args, **kwargs) # numba: DELETE elif fprime_included: # numba: DELETE fval, fder = func(guess, *args, **kwargs) elif fprime2 is not None: #numba: DELETE fval = func(guess, *args, **kwargs) #numba: DELETE fder = fprime(guess, *args, **kwargs) #numba: DELETE fcurv = fprime2(guess, *args, **kwargs) #numba: DELETE else: #numba: DELETE fval = func(guess, *args, **kwargs) #numba: DELETE fder = fprime(guess, *args, **kwargs) #numba: DELETE # print(f"Guess={guess} value={fval} derivative={fder}") # if fprime2_included: # numba: uncomment # fval, fder, fder2 = func(guess, *args) # numba: uncomment # else: # numba: uncomment # fval, fder = func(guess, *args) # numba: uncomment # fder2 = 0.0 # numba: uncomment if require_bracket and fval*fval1 >= 0.0: continue if isnan(fval) or isinf(fval) or fder == 0.0: continue else: p0 = guess recovered = True # print(f'Found point the search can continue from {p0}') break if not recovered: raise UnconvergedError("Derivative became zero on iteration %d, guess=%f " %(it, 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*fcurv*fder_inv) p = damping_func(p0, -step, damping) elif fprime2 is None or isnan(fcurv) or isinf(fcurv): p = p0 - step*damping else: newton_step = step*damping halley_step = newton_step/(1.0 - 0.5*step*fcurv*fder_inv) if halley_step*newton_step < 0.0: # Both changes go in a different direction, therefore only use the newton step p = p0 - newton_step else: p = p0 - halley_step if bisection and a is not None and b is not None: if ((not (a < p < b) and not (b < p < a)) or it > 40): # Arbitrary switch to bisection if we should have converged by now. # A common issue is when the derivative has a small numerical issue around something # if p < 0.0: # if p < a: # print('bisecting') p = 0.5*(a + b) # else: # if p > b: # 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 ytol is not None and 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 ytol is not None and 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 (unknown fval) it += 1 if ytol is not None and xtol is not None: # Meet both tolerance - new value is under ytol, and old value ytol_met = abs(fval) <= ytol xtol_met = abs(p - p0) <= abs(xtol*p) # Has to be less or equal to for the case of zero! if (xtol_met and ytol_met) or (not require_xtol and ytol_met): 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, fval2, fval3 = fval, fval0, fval1, fval2 if isnan(p) or isinf(p): raise UnconvergedError("Cannot continue - math error in function value on iteration %d, guess=%f, value=%f"%(it, p, fval)) # print([p, p0, p1, p2], [fval, fval0, fval1, fval2], ytol_met, xtol_met) p0, p1, p2, p3 = p, p0, p1, p2 raise UnconvergedError("Failed to converge; maxiter (%d) reached, point=%g, error=%g" %(maxiter, p, fval)) 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 check_jacobian(x, j, func, jac_error_allowed=False): '''Basic function to check a jacobian for inf and nan, and approximate it if it does contain those elements. ''' bad_jacobian = False for jrow in j: for jval in jrow: if isinf(jval) or isnan(jval): bad_jacobian = True break if bad_jacobian: break if jac_error_allowed: j = jacobian(func, x, scalar=False) bad_jacobian = False for jrow in j: for jval in jrow: if isinf(jval) or isnan(jval): bad_jacobian = True break if bad_jacobian: break if bad_jacobian: raise ValueError("Cannot continue - nan or inf error found in Jacobian") return j def newton_system(f, x0, jac, xtol=None, ytol=None, maxiter=100, damping=1.0, args=(), damping_func=None, line_search=False, require_progress=False, check_numbers=False, Armijo=False, Armijo_c1=1e-4, solve_func=py_solve, with_point=False, jac_error_allowed=False): # numba: delete # solve_func=np.linalg.solve): # 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 N = len(fcur) 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 x = x0 if not jac_also: # numba: delete j = jac(x, *args) # numba: delete if check_numbers:# numba: delete j = check_jacobian(x=x0, j=j, func=f, jac_error_allowed=jac_error_allowed)# numba: delete factors = newton_line_search_factors if line_search else newton_line_search_factors_disabled jac_updated = True # numba: delete iteration = 1 while iteration < maxiter: try:# numba: delete dx = solve_func(j, [-v for v in fcur]) # numba: delete except Exception as e:# numba: delete if jac_error_allowed:# numba: delete j = jacobian(f, x, scalar=False) # numba: delete dx = solve_func(j, [-v for v in fcur]) # numba: delete else:# numba: delete raise e# numba: delete # dx = solve_func(j, -fcur) # numba: uncomment for factor in factors: # print(factor) 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: # numba: delete if jac_also: # numba: delete fnew, jnew = f(xnew, *args) # numba: delete jac_updated = True # numba: delete else: # numba: delete fnew = f(xnew, *args) # numba: delete # print(xnew, 'xnew') # print(fnew, 'fnew') do_next = False # numba: delete if check_numbers: # numba: delete for v in fnew: # numba: delete if isinf(v) or isnan(v): # numba: delete do_next = True # numba: delete if do_next and factor != factors[-1]:# numba: delete # print('math error, trying next point') continue# numba: delete except: # numba: delete # print(f'Line search calculation with point failed') # numba: delete continue # numba: delete # fnew, jnew = f(xnew, *args) # numba: uncomment err_new = 0.0 for v in fnew: err_new += abs(v) if check_numbers: for v in fnew: if isinf(v) or isnan(v): raise ValueError("Cannot continue - math error in function value") # print(f'Line search with error={err_new}, factor {mult}' ) if err_new < err0: if Armijo: phi = 0.0 for v in fnew: phi += v*v Armijo_lhs = phi phi0 = 0.0 for v in fcur: phi0 += v*v derphi0 = 0.0 for i in range(N): temp = 0.0 for jidx in range(N): temp += fcur[jidx]*j[jidx][i] # Armijo_grad[i] = 2.0*temp derphi0 += 4.0*temp*temp # derphi0 = float(np.dot(Armijo_grad, Armijo_grad)) Armijo_rhs = Armijo_c1*factor*derphi0 + phi0 if Armijo_lhs <= Armijo_rhs: if not jac_also: # numba: delete jac_updated = False # numba: delete break # else: # print('lhs and rhs not same', Armijo_lhs, Armijo_rhs) else: if not jac_also: # numba: delete jac_updated = False # numba: delete break if (line_search and err_new > err0) and require_progress: raise ValueError("Completed line search without reducing the objective function error, cannot proceed") fcur = fnew if err_new >= err0: # edge case, didn't make any progress but keep going if require_progress: raise ValueError("Completed line search without reducing the objective function error, cannot proceed") else:# numba: delete if not jac_also: # numba: delete jnew = jac(xnew, *args) # numba: delete jac_updated = True # numba: delete # print(xnew) 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_updated: # numba: delete jnew = jac(x, *args) # numba: delete if check_numbers:# numba: delete jnew = check_jacobian(x=xnew, j=jnew, func=f, jac_error_allowed=jac_error_allowed) # numba: delete # print('jnew', jnew) j = jnew jac_updated = False # 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") if with_point: # numba: delete return x, iteration, fcur, j # numba: delete 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 homotopy_function(x, lambd, F0, objf, args=()): ''' .. math:: H(x, lambda) = F(x) - (1-lambda)*F(x_original) x : list[float] The current guesses, [-] lambd : float Lambda parameter, between 0 and 1, [-] F0 : list[float] Original error functions from the very first iteration, [-] objf: callable Error function, [-] args : tuple Extra functions to objevtive function, [-] ''' Fx = objf(x, *args) const = 1.0 - lambd # print(f'x={x}, Fx={Fx}, F0={F0}, err={[Fx[i] - const*F0[i] for i in range(len(F0))]}') return [Fx[i] - const*F0[i] for i in range(len(F0))] def homotopy_solver(f, x0, jac, xtol=1e-8, ytol=None, maxiter=100, lambd_step=0.1, damping=1.0, args=(), line_search=True, terminal_xtol=1e-10, terminal_ytol=None, terminal_lambd_criteria=0.9, solve_func=py_solve): f0 = f(x0, *args) working_lambda_step = lambd_step lambd_previous = lambd = 0.0 N = len(f0) remaining_iters = maxiter homotopy_iter = 0 failed_iters = 0 f_n = lambda x, *discarded_args: homotopy_function(x, lambd, f0, f, args=args) while abs(lambd - 1) > 1e-15: lambd += working_lambda_step if lambd > 1.0: lambd = 1.0 # Guess a next step by computing dx_dlambda with a numerical # evaluation of the derivative with respect to lambda # using the previous lambda value at the current point delta_lambd = working_lambda_step if homotopy_iter == 0: Fx = f0 J = jac(x0, *args) H0 = [Fx[i] - (1.0 - lambd_previous)*f0[i] for i in range(N)] H1 = [Fx[i] - (1.0 - lambd)*f0[i] for i in range(N)] const = -1.0/delta_lambd dx_dlambd = solve_func(J, [const*(H1[i] - H0[i]) for i in range(N)]) # This gets the x values a little closer to the new H function x1 = [x0[i] + delta_lambd*dx_dlambd[i] for i in range(N)] # Now solve the objective function using a newton solver with line search use_xtol, use_ytol = (xtol, ytol) if lambd < terminal_lambd_criteria else (terminal_xtol, terminal_ytol) try: # last objf and J call is wasted- would be nice to reuse them in the earlier code x2, niter, Fx, J = newton_system(f=f_n, x0=x1, jac=jac, xtol=use_xtol, ytol=use_ytol, maxiter=remaining_iters, damping=damping, args=args, solve_func=solve_func, line_search=line_search, with_point=True) except: # undo the lambda change, continue the loop with half the lambda step lambd -= working_lambda_step working_lambda_step *= 0.5 failed_iters += 1 if failed_iters > 5: raise UnconvergedError("Failed to converge") continue remaining_iters -= niter delta_lambd = lambd_step x0 = x2 lambd_previous = lambd # print(f'Converged lambda step {lambd}') homotopy_iter += 1 return x0, maxiter - remaining_iters 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 = matrix_vector_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 = matrix_vector_dot(J, z) d = [-i for i in s] dmu = [d[i]-u[i] for i in eqns] dmu_d = dot_product(dmu, d) den_inv = 1.0/dot_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 fixed_point(f, x0, xtol=None, ytol=None, maxiter=100, damping=1.0, args=(), require_progress=False, check_numbers=False): fcur = f(x0, *args) 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.0 x = [xi - fi for xi, fi in zip(x0, fcur)] iteration = 1 while iteration < maxiter: fcur = f(x, *args) if check_numbers: for v in fcur: if isinf(v) or isnan(v): raise ValueError("Cannot continue - math error in function value") x = [xi - fi*damping for xi, fi in zip(x, fcur)] err1 = 0.0 for v in fcur: err1 += abs(v) iteration += 1 if xtol is not None: if (norm2(fcur) < xtol) and (ytol is None or err1 < ytol): break elif ytol is not None: if err1 < ytol: break if err1 >= err0 and require_progress: raise ValueError("Fixed point is not making progress, cannot proceed") err0 = err1 if xtol is not None and norm2(fcur) > xtol: raise UnconvergedError("Failed to converge") if ytol is not None: err0 = 0.0 for v in fcur: err0 += abs(v) if err0 > ytol: raise UnconvergedError("Failed to converge") return x, iteration def aitken_delta_squared_accelerate(x, x1, x2, max_change_ratio=1.0): N = len(x) out = [0.0]*N for i in range(N): dx1 = x2[i] - x1[i] dx = x1[i] - x[i] den = dx1 - dx if den != 0.0: change = -dx1*dx1/den # print(i, abs(change/x[i])) if change != 0.0 and x[i] != 0.0 and abs(change/x[i]) < max_change_ratio: out[i] = change return out def fixed_point_aitken(f, x0, xtol=None, ytol=None, maxiter=100, damping=1.0, args=(), require_progress=False, check_numbers=False, acc_frequency=4, acc_damping=1.0, acc_max_change_ratio=1.0): all_guesses = [x0] fcur = f(x0, *args) 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.0 x = [xi - fi for xi, fi in zip(x0, fcur)] iteration = 1 while iteration < maxiter: all_guesses.append(x) fcur = f(x, *args) if check_numbers: for v in fcur: if isinf(v) or isnan(v): raise ValueError("Cannot continue - math error in function value") if (iteration % acc_frequency) == 0: x1, x2 = all_guesses[-2], all_guesses[-3] dx = aitken_delta_squared_accelerate(x, x1, x2, acc_max_change_ratio) x = [xi + dxi*acc_damping for xi, dxi in zip(x, dx)] else: x = [xi - fi*damping for xi, fi in zip(x, fcur)] err1 = 0.0 for v in fcur: err1 += abs(v) iteration += 1 if xtol is not None: if (norm2(fcur) < xtol) and (ytol is None or err1 < ytol): break elif ytol is not None: if err1 < ytol: break if err1 >= err0 and require_progress: raise ValueError("Fixed point is not making progress, cannot proceed") err0 = err1 if xtol is not None and norm2(fcur) > xtol: raise UnconvergedError("Failed to converge") if ytol is not None: err0 = 0.0 for v in fcur: err0 += abs(v) if err0 > ytol: raise UnconvergedError("Failed to converge") return x, iteration def dem(x, x1, x2, max_change_ratio=1.0): N = len(x) dx1 = [x1[i] - x2[i] for i in range(N)] dx = [x[i] - x1[i] for i in range(N)] dot_dx1_dx = 0.0 dot_dx_dx = 0.0 for i in range(N): dot_dx_dx += dx[i]*dx[i] dot_dx1_dx += dx1[i]*dx[i] if dot_dx_dx == 0: return [0.0]*N lbda = dot_dx1_dx/dot_dx_dx den = 1.0 - lbda if den == 0.0: return [0.0]*N # Not entirely sure why negative is needed den_inv = abs(1.0/den) dacc = [dx[i]*den_inv for i in range(N)] return dacc def gdem(x, x1, x2, x3, max_change_ratio=1.0): N = len(x) dx2 = [x[i] - x3[i] for i in range(N)] dx1 = [x[i] - x2[i] for i in range(N)] dx = [x[i] - x1[i] for i in range(N)] b01, b02, b12, b11, b22 = 0.0, 0.0, 0.0, 0.0, 0.0 for i in range(N): b01 += dx[i]*dx1[i] b02 += dx[i]*dx2[i] b12 += dx1[i]*dx2[i] b11 += dx1[i]*dx1[i] b22 += dx2[i]*dx2[i] to_den = b11*b22 - b12*b12 if abs(to_den) == 0: return dem(x, x1, x2, max_change_ratio) den_inv = 1.0/(to_den) mu1 = den_inv*(b02*b12 - b01*b22) mu2 = den_inv*(b01*b12 - b02*b11) to_den2 = 1.0 + mu1 + mu2 if to_den2 == 0: return dem(x, x1, x2, max_change_ratio) factor = 1.0/(to_den2) dx_out = [dxi for dxi in dx] for i in range(N): change = factor*(dx[i] - mu2*dx1[i]) if x[i] != 0.0 and abs(change/x[i]) < max_change_ratio: dx_out[i] = change return dx_out def fixed_point_gdem(f, x0, xtol=None, ytol=None, maxiter=100, damping=1.0, args=(), require_progress=False, check_numbers=False, acc_frequency=5, acc_damping=1.0, acc_max_change_ratio=1.0): all_guesses = [x0] fcur = f(x0, *args) 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.0 x = [xi - fi for xi, fi in zip(x0, fcur)] iteration = 1 while iteration < maxiter: all_guesses.append(x) fcur = f(x, *args) if check_numbers: for v in fcur: if isinf(v) or isnan(v): raise ValueError("Cannot continue - math error in function value") if (iteration % acc_frequency) == 0: x1, x2, x3 = all_guesses[-2], all_guesses[-3], all_guesses[-4] dx = gdem(x, x1, x2, x3, acc_max_change_ratio) # print(np.array([xi - fi*acc_damping for xi, fi in zip(x, fcur)])/[xi + dxi*acc_damping for xi, dxi in zip(x, dx)]) x = [xi + dxi*acc_damping for xi, dxi in zip(x, dx)] else: x = [xi - fi*damping for xi, fi in zip(x, fcur)] err1 = 0.0 for v in fcur: err1 += abs(v) iteration += 1 if xtol is not None: if (norm2(fcur) < xtol) and (ytol is None or err1 < ytol): break elif ytol is not None: if err1 < ytol: break if err1 >= err0 and require_progress: raise ValueError("Fixed point is not making progress, cannot proceed") err0 = err1 if xtol is not None and norm2(fcur) > xtol: raise UnconvergedError("Failed to converge") if ytol is not None: err0 = 0.0 for v in fcur: err0 += abs(v) if err0 > ytol: raise UnconvergedError("Failed to converge") return x, iteration # def compute_accelerated_step( # residuals: List[List[float]], # x_hist: List[List[float]], # gx_hist: List[List[float]], # reg: float = 1e-8, # mixing_param: float = 1.0 # ) -> List[float]: # """Compute Anderson acceleration coefficients and the accelerated iterate.""" # # Compute R = Ft @ Ft.T with regularization # RR = matrix_multiply(residuals, transpose(residuals)) # N = len(RR) # for i in range(N): # RR[i][i] += reg # try: # RR_inv = inv(RR) # alpha = sum_matrix_rows(RR_inv) # except: # # Fallback to least squares if matrix is singular # ones = [1.0] * len(residuals) # alpha = gelsd(RR, ones)[0] # # Normalize alpha # alpha_sum = sum(alpha) # # DO NOT REMOVE part of what is needed to switch to f in residual_to_fixed_point form # if alpha_sum == 0: # raise ValueError("Sum of alpha coefficients is zero") # elif alpha_sum < 0: # alpha = [-a for a in alpha] # alpha_sum = -alpha_sum # # Compute the accelerated iterate # dim = len(x_hist[0]) # x_acc = [0.0] * dim # for a, x, gx in zip(alpha, x_hist, gx_hist): # for i in range(dim): # x_acc[i] += (1 - mixing_param) * a * x[i] + mixing_param * a * gx[i] # return x_acc # def anderson_step( # x_hist: List[List[float]], # gx_hist: List[List[float]], # residuals_hist: List[List[float]], # x: List[float], # window_size: int, # reg: float = 1e-8, # mixing_param: float = 1.0 # ) -> Tuple[List[float], List[List[float]], List[List[float]], List[List[float]]]: # """Perform one step of Anderson acceleration.""" # # First iteration case # if not x_hist: # return x, [x], [], [] # # Update histories # x_prev = x_hist[-1] # residual = [xi - xp for xi, xp in zip(x, x_prev)] # new_residuals_hist = residuals_hist + [residual] # new_gx_hist = gx_hist + [x] # # Maintain window size # if len(new_residuals_hist) > window_size + 1: # new_residuals_hist = new_residuals_hist[1:] # new_gx_hist = new_gx_hist[1:] # # Compute accelerated iterate # x_acc = compute_accelerated_step( # new_residuals_hist, # x_hist, # new_gx_hist, # reg, # mixing_param # ) # # Update x history # new_x_hist = x_hist + [x_acc] # if len(new_x_hist) > window_size + 1: # new_x_hist = new_x_hist[1:] # return x_acc, new_x_hist, new_gx_hist, new_residuals_hist # def fixed_point_anderson( # f: Callable, # x0: List[float], # xtol: float = 1e-7, # ytol: float = None, # maxiter: int = 100, # args: tuple = (), # require_progress: bool = False, # check_numbers: bool = False, # window_size: int = 5, # reg: float = 1e-8, # mixing_param: float = 1.0 # ) -> Tuple[List[float], int]: # """Fixed point iteration with Anderson acceleration (functional version).""" # # Initialize state # x_hist: List[List[float]] = [] # gx_hist: List[List[float]] = [] # residuals_hist: List[List[float]] = [] # x = x0 # fcur = f(x, *args) # # Check initial convergence for ytol # err0 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 # if ytol is not None and xtol is None and err0 < ytol: # return x0, 0 # # Main iteration loop # for iteration in range(maxiter): # x_new = f(x, *args) # # Check for inf/nan # if check_numbers and any(isnan(v) or isinf(v) for v in x_new): # raise ValueError("Cannot continue - math error in function value") # # Apply Anderson acceleration # x_acc, x_hist, gx_hist, residuals_hist = anderson_step( # x_hist, gx_hist, residuals_hist, x_new, # window_size, reg, mixing_param # ) # # Calculate errors # err1 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 # # Check progress # if require_progress and ytol is not None and err1 >= err0: # raise ValueError("Fixed point is not making progress") # # Update error # err0 = err1 if ytol is not None else 0.0 # # Check convergence # if xtol is not None: # x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) # if x_err < xtol and (ytol is None or err1 < ytol): # return x_acc, iteration # elif ytol is not None and err1 < ytol: # return x_acc, iteration # x = x_acc # # Check final convergence # x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) # if xtol is not None and x_err > xtol: # raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {x_err}") # if ytol is not None and err1 > ytol: # raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {err1}") # return x, iteration def compute_accelerated_step( residuals, x_hist, gx_hist, reg: float = 1e-8, mixing_param: float = 1.0 ): """Compute Anderson acceleration coefficients and the accelerated iterate.""" # Compute R = Ft @ Ft.T with regularization RR = matrix_multiply(residuals, transpose(residuals)) N = len(RR) for i in range(N): RR[i][i] += reg try: RR_inv = inv(RR) alpha = sum_matrix_rows(RR_inv) except: # Fallback to least squares if matrix is singular ones = [1.0] * len(residuals) alpha = gelsd(RR, ones)[0] # Normalize alpha alpha_sum = sum(alpha) # DO NOT REMOVE part of what is needed to switch to f in residual_to_fixed_point form if alpha_sum == 0: raise ValueError("Sum of alpha coefficients is zero") elif alpha_sum < 0: alpha = [-a for a in alpha] alpha_sum = -alpha_sum # sometimes alpha needs to be negative alpha = [a / alpha_sum for a in alpha] # Compute the accelerated iterate dim = len(x_hist[0]) x_acc = [0.0] * dim for a, x, gx in zip(alpha, x_hist, gx_hist): for i in range(dim): # Flip between these to change the basis # x_acc[i] += (1 - mixing_param) * a * x[i] + mixing_param * a * gx[i] x_acc[i] +=(1 - mixing_param) * a * x[i] - mixing_param * a * gx[i] return x_acc # def anderson_step( # x_hist: List[List[float]], # gx_hist: List[List[float]], # residuals_hist: List[List[float]], # x: List[float], # window_size: int, # reg: float = 1e-8, # mixing_param: float = 1.0 # ) -> Tuple[List[float], List[List[float]], List[List[float]], List[List[float]]]: # """Perform one step of Anderson acceleration.""" # # First iteration case # if not x_hist: # return x, [x], [], [] # # Update histories # x_prev = x_hist[-1] # # Flip between these to change the basis # residual = [xi - xp for xi, xp in zip(x, x_prev)] # # residual = [xi + xp for xi, xp in zip(x, x_prev)] # new_residuals_hist = residuals_hist + [residual] # new_gx_hist = gx_hist + [x] # # Maintain window size # if len(new_residuals_hist) > window_size + 1: # new_residuals_hist = new_residuals_hist[1:] # new_gx_hist = new_gx_hist[1:] # # Compute accelerated iterate # x_acc = compute_accelerated_step( # new_residuals_hist, # x_hist, # new_gx_hist, # reg, # mixing_param # ) # # Update x history # new_x_hist = x_hist + [x_acc] # if len(new_x_hist) > window_size + 1: # new_x_hist = new_x_hist[1:] # return x_acc, new_x_hist, new_gx_hist, new_residuals_hist # def fixed_point_anderson( # f: Callable, # x0: List[float], # xtol: float = 1e-7, # ytol: float = None, # maxiter: int = 100, # args: tuple = (), # require_progress: bool = False, # check_numbers: bool = False, # window_size: int = 5, # reg: float = 1e-8, # mixing_param: float = 1.0 # ) -> Tuple[List[float], int]: # """Fixed point iteration with Anderson acceleration (functional version).""" # # Initialize state # x_hist: List[List[float]] = [] # gx_hist: List[List[float]] = [] # residuals_hist: List[List[float]] = [] # x = x0 # # f = residual_to_fixed_point(f) # this commented out means the first lines should be uncommented; second lines commented # fcur = f(x, *args) # # Check initial convergence for ytol # err0 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 # if ytol is not None and xtol is None and err0 < ytol: # return x0, 0 # # Main iteration loop # for iteration in range(maxiter): # x_new = f(x, *args) # # Check for inf/nan # if check_numbers and any(isnan(v) or isinf(v) for v in x_new): # raise ValueError("Cannot continue - math error in function value") # # Apply Anderson acceleration # x_acc, x_hist, gx_hist, residuals_hist = anderson_step( # x_hist, gx_hist, residuals_hist, x_new, # window_size, reg, mixing_param # ) # # Calculate errors # err1 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 # # Check progress # if require_progress and ytol is not None and err1 >= err0: # raise ValueError("Fixed point is not making progress") # # Update error # err0 = err1 if ytol is not None else 0.0 # # Check convergence # if xtol is not None: # x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) # if x_err < xtol and (ytol is None or err1 < ytol): # return x_acc, iteration # elif ytol is not None and err1 < ytol: # return x_acc, iteration # x = x_acc # # Check final convergence # x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) # if xtol is not None and x_err > xtol: # raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {x_err}") # if ytol is not None and err1 > ytol: # raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {err1}") # return x, iteration # def anderson_step( # x_hist: List[List[float]], # gx_hist: List[List[float]], # residuals_hist: List[List[float]], # x: List[float], # window_size: int, # reg: float = 1e-8, # mixing_param: float = 1.0, # current_iteration: int = 0, # delay_iterations: int = 0 # ) -> Tuple[List[float], List[List[float]], List[List[float]], List[List[float]]]: # """Perform one step of Anderson acceleration with history building during delay.""" # # First iteration case # if not x_hist: # return x, [x], [x], [] # # Update histories # x_prev = x_hist[-1] # # Flip between these to change the basis # residual = [xi - xp for xi, xp in zip(x, x_prev)] # # residual = [xi + xp for xi, xp in zip(x, x_prev)] # new_x_hist = x_hist + [x] # new_gx_hist = gx_hist + [x] # new_residuals_hist = residuals_hist + [residual] # # Maintain window size # if len(new_residuals_hist) > window_size + 1: # new_residuals_hist = new_residuals_hist[1:] # new_gx_hist = new_gx_hist[1:] # new_x_hist = new_x_hist[1:] # # During delay period, just return the current point but keep updating histories # if current_iteration < delay_iterations: # return x, new_x_hist, new_gx_hist, new_residuals_hist # # After delay, compute accelerated iterate using accumulated history # x_acc = compute_accelerated_step( # new_residuals_hist, # new_x_hist, # new_gx_hist, # reg, # mixing_param # ) # # Update x history with accelerated point # final_x_hist = new_x_hist[:-1] + [x_acc] # return x_acc, final_x_hist, new_gx_hist, new_residuals_hist # def fixed_point_anderson( # f: Callable, # x0: List[float], # xtol: float = 1e-7, # ytol: float = None, # maxiter: int = 100, # args: tuple = (), # require_progress: bool = False, # check_numbers: bool = False, # window_size: int = 5, # reg: float = 1e-8, # mixing_param: float = 1.0, # delay_iterations: int = 0 # ) -> Tuple[List[float], int]: # """Fixed point iteration with Anderson acceleration and history building during delay.""" # # Initialize state # x_hist: List[List[float]] = [] # gx_hist: List[List[float]] = [] # residuals_hist: List[List[float]] = [] # x = x0 # # f = residual_to_fixed_point(f) # this commented out means the first lines should be uncommented; second lines commented # fcur = f(x, *args) # # Check initial convergence for ytol # err0 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 # if ytol is not None and xtol is None and err0 < ytol: # return x0, 0 # # Main iteration loop # for iteration in range(maxiter): # x_new = f(x, *args) # # Check for inf/nan # if check_numbers and any(isnan(v) or isinf(v) for v in x_new): # raise ValueError("Cannot continue - math error in function value") # # Apply Anderson acceleration with delay but keep building history # x_acc, x_hist, gx_hist, residuals_hist = anderson_step( # x_hist, gx_hist, residuals_hist, x_new, # window_size, reg, mixing_param, # iteration, delay_iterations # ) # # Calculate errors # err1 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 # # Check progress # if require_progress and ytol is not None and err1 >= err0: # raise ValueError("Fixed point is not making progress") # # Update error # err0 = err1 if ytol is not None else 0.0 # # Check convergence # if xtol is not None: # x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) # if x_err < xtol and (ytol is None or err1 < ytol): # return x_acc, iteration # elif ytol is not None and err1 < ytol: # return x_acc, iteration # x = x_acc # # Check final convergence # x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) # if xtol is not None and x_err > xtol: # raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {x_err}") # if ytol is not None and err1 > ytol: # raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {err1}") # return x, iteration def compute_accelerated_step( residuals, x_hist, gx_hist, reg: float = 1e-8, mixing_param: float = 1.0 ): """Compute Anderson acceleration coefficients and the accelerated iterate.""" # Compute R = Ft @ Ft.T with regularization RR = matrix_multiply(residuals, transpose(residuals)) N = len(RR) for i in range(N): RR[i][i] += reg try: RR_inv = inv(RR) alpha = sum_matrix_rows(RR_inv) except: # Fallback to least squares if matrix is singular ones = [1.0] * len(residuals) alpha = gelsd(RR, ones)[0] # Normalize alpha alpha_sum = sum(alpha) if alpha_sum == 0: raise ValueError("Sum of alpha coefficients is zero") elif alpha_sum < 0: alpha = [-a for a in alpha] alpha_sum = -alpha_sum # sometimes alpha needs to be negative alpha = [a / alpha_sum for a in alpha] # Compute the accelerated iterate dim = len(x_hist[0]) x_acc = [0.0] * dim for a, x, gx in zip(alpha, x_hist, gx_hist): for i in range(dim): # Flip between these to change the basis x_acc[i] += (1 - mixing_param) * a * x[i] + mixing_param * a * gx[i] # x_acc[i] +=(1 - mixing_param) * a * x[i] - mixing_param * a * gx[i] return x_acc def anderson_step( x_hist, gx_hist, residuals_hist, x, window_size: int, reg: float = 1e-8, mixing_param: float = 1.0 ): """Perform one step of Anderson acceleration.""" # First iteration case if not x_hist: return x, [x], [], [] # Update histories x_prev = x_hist[-1] # Flip between these to change the basis residual = [xi - xp for xi, xp in zip(x, x_prev)] # residual = [xi + xp for xi, xp in zip(x, x_prev)] new_residuals_hist = residuals_hist + [residual] new_gx_hist = gx_hist + [x] # Maintain window size if len(new_residuals_hist) > window_size + 1: new_residuals_hist = new_residuals_hist[1:] new_gx_hist = new_gx_hist[1:] # Compute accelerated iterate x_acc = compute_accelerated_step( new_residuals_hist, x_hist, new_gx_hist, reg, mixing_param ) # Update x history new_x_hist = x_hist + [x_acc] if len(new_x_hist) > window_size + 1: new_x_hist = new_x_hist[1:] return x_acc, new_x_hist, new_gx_hist, new_residuals_hist def fixed_point_anderson( f, x0, xtol: float = 1e-7, ytol: float = None, maxiter: int = 100, args: tuple = (), require_progress: bool = False, check_numbers: bool = False, window_size: int = 5, reg: float = 1e-8, mixing_param: float = 1.0, initial_iterations: int = 4, damping: float = 0.9, acc_damping: float = 1, # Damping for acceleration phase max_step_size: float = 10000, # Maximum allowed relative step size phase_in_steps: int = 3 # Number of steps to phase in acceleration ): # Initialize state x_hist = [] gx_hist = [] residuals_hist = [] x = x0 # f = residual_to_fixed_point(f) # this commented out means the first lines should be uncommented; second lines commented fcur = f(x, *args) # Truncate window size based on the dimensionality of the problem # sometimes being 1 larger can be OK, for now not allowing it if len(fcur) < window_size: window_size = len(fcur) # Check initial convergence for ytol err0 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 if ytol is not None and xtol is None and err0 < ytol: return x0, 0 # Main iteration loop for iteration in range(maxiter): x_new = f(x, *args) # Check for inf/nan if check_numbers and any(isnan(v) or isinf(v) for v in x_new): raise ValueError("Cannot continue - math error in function value") if iteration < initial_iterations: # Damped fixed-point iteration x_acc = [ (1 - damping) * x[i] + damping * x_new[i] for i in range(len(x)) ] # Collect history for later use if iteration > 0: residual = [xi - xp for xi, xp in zip(x_new, x)] residuals_hist.append(residual) gx_hist.append(x_new) x_hist.append(x_acc) # Maintain window size if len(residuals_hist) > window_size: residuals_hist = residuals_hist[1:] gx_hist = gx_hist[1:] x_hist = x_hist[1:] else: # Apply Anderson acceleration with safeguards x_acc_raw, x_hist, gx_hist, residuals_hist = anderson_step( x_hist, gx_hist, residuals_hist, x_new, window_size, reg, mixing_param ) # Calculate phase-in factor (gradually increase from 0 to 1) phase = min(1.0, (iteration - initial_iterations + 1) / phase_in_steps) # Compute damped acceleration step x_acc_damped = [] for i in range(len(x)): # Regular damped fixed-point step fp_step = (1 - damping) * x[i] + damping * x_new[i] # Anderson acceleration step with its own damping acc_step = (1 - acc_damping) * x[i] + acc_damping * x_acc_raw[i] # Blend between fixed-point and acceleration based on phase x_acc_i = (1 - phase) * fp_step + phase * acc_step # Limit maximum step size relative to current position max_change = abs(x[i] * (1.0-max_step_size)) change = x_acc_i - x[i] if abs(change) > max_change: # Clamp the change to the maximum allowed x_acc_i = x[i] + max_change * (1 if change > 0 else -1) x_acc_damped.append(x_acc_i) x_acc = x_acc_damped x_hist[-1] = x_acc # Calculate errors err1 = sum(abs(v) for v in fcur) if ytol is not None else 0.0 # Check progress if require_progress and ytol is not None and err1 >= err0: raise ValueError("Fixed point is not making progress") # Update error err0 = err1 if ytol is not None else 0.0 # Check convergence if xtol is not None: x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) if x_err < xtol and (ytol is None or err1 < ytol): return x_acc, iteration elif ytol is not None and err1 < ytol: return x_acc, iteration x = x_acc # Check final convergence x_err = max(abs((a - b) / abs(b)) for a, b in zip(x_acc, x)) if xtol is not None and x_err > xtol: raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {x_err}") if ytol is not None and err1 > ytol: raise ValueError(f"Failed to converge after {maxiter} iterations. Error: {err1}") return x, iteration def fixed_point_anderson_residual( f, x0, xtol: float = 1e-7, ytol: float = None, maxiter: int = 100, args: tuple = (), require_progress: bool = False, check_numbers: bool = False, window_size: int = 5, reg: float = 1e-8, initial_iterations: int = 1000, damping: float = 0.3, acc_damping: float = 1 ): x_hist = [] gx_hist = [] # Will store the updates/steps residuals_hist = [] x = x0 # Evaluate initial residual fcur = f(x, *args) if len(fcur) < window_size: window_size = len(fcur) err0 = sum(abs(v) for v in fcur) if ytol is not None and xtol is None and err0 < ytol: return x0, 0 for iteration in range(maxiter): # Compute residual fcur = f(x, *args) # Check for inf/nan if check_numbers and any(isnan(v) or isinf(v) for v in fcur): raise ValueError("Cannot continue - math error in function value") if iteration < initial_iterations: # Damped step during initial phase x_step = [-fi * damping for fi in fcur] x_new = [xi + si for xi, si in zip(x, x_step)] # Collect history (store steps instead of points) if iteration > 0: # For history, store the full steps full_step = [-fi for fi in fcur] residuals_hist.append(fcur) gx_hist.append(full_step) x_hist.append(x) # Maintain window size if len(residuals_hist) > window_size: residuals_hist = residuals_hist[1:] gx_hist = gx_hist[1:] x_hist = x_hist[1:] else: # Full steps for acceleration phase x_step = [-fi for fi in fcur] # Update histories residuals_hist.append(fcur) gx_hist.append(x_step) x_hist.append(x) if len(residuals_hist) > window_size: residuals_hist = residuals_hist[1:] gx_hist = gx_hist[1:] x_hist = x_hist[1:] # Compute Anderson coefficients RR = matrix_multiply(residuals_hist, transpose(residuals_hist)) N = len(RR) for i in range(N): RR[i][i] += reg try: RR_inv = inv(RR) alpha = sum_matrix_rows(RR_inv) except: ones = [1.0] * len(residuals_hist) alpha = gelsd(RR, ones)[0] # Normalize alpha alpha_sum = sum(alpha) if alpha_sum == 0: raise ValueError("Sum of alpha coefficients is zero") alpha = [a / alpha_sum for a in alpha] # Apply acceleration to the steps acc_step = [0.0] * len(x) for a, step in zip(alpha, gx_hist): for i in range(len(x)): acc_step[i] += a * step[i] # Apply accelerated step with damping x_new = [xi + si * acc_damping for xi, si in zip(x, acc_step)] # Calculate errors err1 = sum(abs(v) for v in fcur) # Check progress if require_progress and err1 >= err0: raise ValueError("Not making progress") err0 = err1 # Check convergence if xtol is not None: if err1 < xtol and (ytol is None or err1 < ytol): return x_new, iteration elif ytol is not None and err1 < ytol: return x_new, iteration x = x_new if xtol is not None and err1 > xtol: raise ValueError(f"Failed to converge after {maxiter} iterations") if ytol is not None and err1 > ytol: raise ValueError(f"Failed to converge after {maxiter} iterations") return x, iteration 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 Sequence 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 # if isinstance(x, (float, int, complex)): # x = [x] n = len(t) 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 arg = x 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: return 0.0 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]* (k + 1) hh = [0.0]*k 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 return sp 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 = [[bispev_inner(tx, ty, c, kx, ky, xi, yi) 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): w = [0.0]*(k+1) tb = t[k] n = len(t) h = [0]*6 hh = [0]*5 te = t[n - k - 1] l1 = k + 1 l2 = l1 + 1 arg = x 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 = l1 - k - 1 for j in range(k + 1): w[j] = h[j] return lx, w def bispev_inner(tx, ty, c, kx, ky, x, y): kx1 = kx + 1 ky1 = ky + 1 nky1 = len(ty) - ky1 lx, wx = init_w(tx, kx, x) ly, wy = init_w(ty, ky, y) sp = 0.0 err = 0.0 for i1 in range(kx1): for j1 in range(ky1): l2 = lx*nky1 + ly + i1*nky1 + j1 a = c[l2]*wx[i1]*wy[j1] - err tmp = sp + a err = (tmp - sp) - a sp = tmp return sp 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_val=1000.0): max_err = 0.0 N = len(data) for i in range(N): if data[i] == 0.0: if calc[i] == 0.0: diff = 0.0 else: diff = max_val else: diff = abs((data[i] - calc[i])/data[i]) if diff > max_val: diff = max_val 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): # if points is not None: # numba: delete # if b < a: # numba: delete # a, b = b, a # numba: delete # points = list(sorted(set([p for p in points if a < p < b]))) # numba: delete # points = [a] + points + [b] # numba: delete # area_accumulator = 0.0 # numba: delete # error_accumulator = 0.0 # numba: delete # for i in range(len(points) - 1): # numba: delete # local_a, local_b = points[i], points[i + 1] # numba: delete # local_area, local_error = quad_adaptive(f, local_a, local_b, args, kronrod_points, kronrod_weights, legendre_weights, epsrel, epsabs) # numba: delete # area_accumulator += local_area # numba: delete # error_accumulator += local_error # numba: delete # return area_accumulator, error_accumulator # numba: delete 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 def bounds_clip_naive(x, low, high): if low is not None: for i in range(N): if x[i] < low[i]: x[i] = low[i] if high is not None: for i in range(N): if x[i] > high[i]: x[i] = high[i] return x def sort_nelder_mead_points_numba(sim, fsim): ind = np.argsort(fsim) sim = sim[ind, :] fsim = fsim[ind] return sim, fsim def sort_nelder_mead_points_python(sim, fsim): # inplace for speed sim = [x for _, x in sorted(zip(fsim, sim))] fsim = list(sorted(fsim)) return sim, fsim def nelder_mead(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=100, maxfun=None, initial_simplex=None, adaptive=False, initial_perturbation=0.05, initial_zero_perturbation=0.00025, low=None, high=None): N = len(x0) has_bounds = low is not None or high is not None if adaptive: dim = float(N) rho = 1.0 chi = 1.0 + 2.0/dim psi = 0.75 - 1.0/(2.0*dim) sigma = 1.0 - 1.0/dim else: rho = 1.0 chi = 2.0 psi = 0.5 sigma = 0.5 rhop1 = 1.0 + rho rhochip1 = rho*chi + 1.0 rhopsip1 = rho*psi + 1.0 rhochi = rho*chi rhopsi = rho*psi onempsi = 1.0 - psi # sort_fun = sort_nelder_mead_points_numba # numba: uncomment sort_fun = sort_nelder_mead_points_python # numba: delete sim = [x0.copy() for _ in range(N+1)] # numba: delete # sim = np.zeros((N+1, N), dtype=float) # numba: uncomment if initial_simplex is None: for k in range(N): if x0[k]!= 0: sim[k+1][k] = (1.0 + initial_perturbation)*x0[k] else: sim[k+1][k] = initial_zero_perturbation else: sim = [x.copy() for x in initial_simplex] # numba: delete # sim = np.asfarray(initial_simplex) # numba: uncomment # Make a new list to store the points. Can initialize it to anything, irrelevant fsim = [x0 for _ in range(N+1)] # numba: delete # fsim = np.zeros((N+1, ), dtype=float) # numba: uncomment for k in range(N + 1): fsim[k] = func(sim[k]) # print(sim, fsim) sim, fsim = sort_fun(sim, fsim) iterations = 0 one_N = 1.0/N while (iterations < maxiter): # Do the convergence checks errcheck1 = 0 sim0 = sim[0] for i in range(1, N+1): asim = sim[i] for j in range(N): a_xdiff = abs(sim0[j] - asim[j]) if a_xdiff > errcheck1: errcheck1 = a_xdiff current_fdiff = 0.0 current_best = fsim[0] for i in range(1, N+1): a_fdiff = abs(current_best - fsim[i]) if a_fdiff > current_fdiff: current_fdiff = a_fdiff if (errcheck1 <= xtol and current_fdiff <= ftol): break sim_last = sim[-1] xbar = [0.0]*N for j in range(N): xbarval = 0.0 simj = sim[j] for i in range(N): xbar[i] += simj[i] for i in range(N): xbar[i] *= one_N # xr = rhop1*xbar - rho*sim_last # numba: uncomment xr = [rhop1*xbar[i] - rho*sim_last[i] for i in range(N)] # numba: delete if has_bounds: xr = bounds_clip_naive(xr, low, high) fxr = func(xr) if fxr < fsim[0]: # xe = rhochip1*xbar - rhochi*sim_last # numba: uncomment xe = [rhochip1*xbar[i] - rhochi*sim_last[i] for i in range(N)] # numba: delete if has_bounds: xe = bounds_clip_naive(xe, low, high) fxe = func(xe) if fxe < fxr: sim[-1] = xe fsim[-1] = fxe else: sim[-1] = xr fsim[-1] = fxr else: if fxr < fsim[-2]: sim[-1] = xr fsim[-1] = fxr else: # contraction if fxr < fsim[-1]: # xc = rhopsip1*xbar - rhopsi*sim_last # numba: uncomment xc = [rhopsip1*xbar[i] - rhopsi*sim_last[i] for i in range(N)] # numba: delete if has_bounds: xc = bounds_clip_naive(xc, low, high) fxc = func(xc) sim[-1] = xc fsim[-1] = fxc else: # contraction inside # xcc = onempsi*xbar + psi*sim_last # numba: uncomment xcc = [onempsi*xbar[i] + psi*sim_last[i] for i in range(N)] # numba: delete if has_bounds: xcc = bounds_clip_naive(xcc, low, high) fxcc = func(xcc) sim[-1] = xcc fsim[-1] = fxcc sim, fsim = sort_fun(sim, fsim) iterations += 1 x = sim[0] fval = min(fsim) return x, fval, iterations def fixed_point_to_residual(f_fixed_point): """ Transforms a fixed-point iteration function to a residual-based function. Parameters: - f_fixed_point: Function that takes x and returns the difference x - thing Returns: - A function that outputs residuals: thing - x """ def residual_function(x, *args): # Get the original fixed-point differences (x - thing) fp_diff = f_fixed_point(x, *args) # Calculate the residuals as (thing - x) return [-diff for diff in fp_diff] return residual_function def residual_to_fixed_point(f_residual): """ Transforms a residual-based function to a fixed-point iteration function. Parameters: - f_residual: Function that takes x and returns residuals (thing - x) Returns: - A function that outputs differences for fixed-point: x - thing """ def fixed_point_function(x, *args): # Get the residuals (thing - x) res = f_residual(x, *args) # Calculate the fixed-point differences as (x - thing) return [-r for r in res] return fixed_point_function scipy_root_options = ('hybr', 'lm', 'broyden1', 'broyden2', 'anderson', 'linearmixing', 'diagbroyden', 'excitingmixing', 'krylov', 'df-sane') scipy_root_options_set = frozenset(scipy_root_options) python_solvers = ('newton_system', 'newton_system_line_search', 'newton_system_line_search_progress', 'homotopy_solver', 'broyden2_python', 'fixed_point', 'fixed_point_aitken', 'fixed_point_gdem', 'fixed_point_anderson') python_solvers_set = frozenset(python_solvers) scipy_minimize_options = ('Nelder-Mead', 'Powell', 'CG', 'BFGS', 'Newton-CG', 'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP', 'trust-constr', 'trust-krylov',# needs hessian 'trust-ncg',# needs hessian 'dogleg', # needs hessian 'trust-exact', # needs hessian ) scipy_requires_hessian_options = ('trust-krylov', 'trust-ncg', 'dogleg', 'trust-exact') scipy_requires_hessian_options_set = frozenset(scipy_requires_hessian_options) scipy_requires_jacobian_options = ('BFGS', 'Newton-CG', 'dogleg', 'trust-ncg', 'trust-krylov') scipy_requires_jacobian_options_set = frozenset(scipy_requires_jacobian_options) scipy_minimize_options_set = frozenset(scipy_minimize_options) multivariate_solvers = python_solvers + scipy_root_options + scipy_minimize_options jacobian_methods = ('python', 'numdifftools_forward', 'numdifftools_reverse', 'numdifftools_central', 'jacobi_forward', 'jacobi_central', 'jacobi_backward') python_jacobians_set = frozenset(['python']) hessian_methods = ('python',) python_hessians_set = frozenset(hessian_methods) class SolverInterface(object): __slots__ = ('minimizing', 'fval_iter', 'jac_iter', 'jac_fval_count', 'method', 'objf_original', 'original_jac', 'xtol', 'ytol', 'maxiter', 'damping', 'jacobian_perturbation', 'jacobian_zero_offset', 'jacobian_method', 'jacobian_order', 'objf_numpy', 'matrix_solver', 'jacobian_numpy', 'solver_numpy', 'objf', 'solver_analytical_jac', 'hess_iter', 'hess_fval_count', 'hessian_method', 'hessian_numpy') def objf_counting(self, *args): self.fval_iter += 1 return self.objf_original(*args) def objf_python_return_numpy(self, x, *args): # function only knows python, solver knows numpy self.fval_iter += 1 return np.array(self.objf_original(x if type(x) is list else x.tolist(), *args)) def objf_numpy_return_python(self, x, *args): # function only knows numpy, solver knows python self.fval_iter += 1 return self.objf_original(np.array(x), *args).tolist() def objf_numpy_minimizing(self, x, *args): self.fval_iter += 1 errs = self.objf_original(np.array(x), *args) tot = 0.0 for v in errs: tot += v*v return float(tot) def objf_python_minimizing(self, x, *args): self.fval_iter += 1 errs = self.objf_original(x if type(x) is list else x.tolist(), *args) tot = 0.0 for v in errs: tot += v*v return tot def jac_minimizing(self, x, *args): fval = self.objf_original(x, *args) jval = self.original_jac(x, *args) N = len(fval) small_jac = [0.0]*N for i in range(N): temp = 0.0 for j in range(N): temp += fval[j]*jval[j][i] small_jac[i] = 2.0*temp return small_jac def jac_minimizing_numpy(self, x, *args): fval = self.objf_original(x, *args) jval = self.original_jac(x, *args) return 2.0*np.dot(fval, jval) def __init__(self, method, objf, jac=None, xtol=1e-8, ytol=None, maxiter=100, damping=1.0, jacobian_method='python', jacobian_perturbation=1e-9, jacobian_zero_offset=1e-7, hessian_method='python', jacobian_order=1, objf_numpy=False, matrix_solver=py_solve): self.method, self.objf_original, self.original_jac = method, objf, jac self.xtol, self.ytol, self.maxiter, self.damping = xtol, ytol, maxiter, damping (self.jacobian_perturbation, self.jacobian_zero_offset, self.jacobian_method, self.jacobian_order, self.hessian_method) = (jacobian_perturbation, jacobian_zero_offset, jacobian_method, jacobian_order, hessian_method) self.objf_numpy, self.matrix_solver = objf_numpy, matrix_solver if jacobian_method == 'analytical': self.jacobian_numpy = objf_numpy else: self.jacobian_numpy = jacobian_method not in python_jacobians_set if hessian_method == 'analytical': self.hessian_numpy = objf_numpy else: self.hessian_numpy = hessian_method not in python_hessians_set # whether or not the solver uses numpy self.solver_numpy = solver_numpy = method not in python_solvers_set self.minimizing = False self.fval_iter = self.jac_iter = self.jac_fval_count = self.hess_iter = self.hess_fval_count = 0 # whether or not the objf uses numpy # if jac is provided it is assumed it is in the same basis if method in scipy_minimize_options_set: self.minimizing = True if jacobian_method == 'scipy' and method in scipy_requires_jacobian_options_set: self.jacobian_method = 'python' self.objf = self.objf_numpy_minimizing if objf_numpy else self.objf_python_minimizing elif solver_numpy: if not objf_numpy: self.objf = self.objf_python_return_numpy else: self.objf = self.objf_counting else: if objf_numpy: self.objf = self.objf_numpy_return_python else: self.objf = self.objf_counting if self.minimizing: if objf_numpy: self.solver_analytical_jac = self.jac_minimizing_numpy else: self.solver_analytical_jac = self.jac_minimizing else: self.solver_analytical_jac = self.original_jac def hessian(self, x, base=None, args=()): self.hess_iter += 1 fval_iter, hessian_method, objf_numpy, hessian_numpy = self.fval_iter, self.hessian_method, self.objf_numpy, self.hessian_numpy return_numpy = type(x) is not list if hessian_method == 'analytical': if objf_numpy and not return_numpy: x = np.array(x) elif return_numpy and not objf_numpy: x = x.tolist() raise NotImplementedError h = self.solver_analytical_hess(x, *args) else: # if the hessian method doesn't speak numpy, convert x to a list if not hessian_numpy: x = x if type(x) is list else x.tolist() else: x = np.array(x) # If the objf doesn't speak numpy but the hessian does, use the converter if self.minimizing: if objf_numpy: objf = self.objf_numpy_minimizing else: objf = self.objf_python_minimizing else: if not hessian_numpy: if objf_numpy: objf = self.objf_numpy_return_python else: objf = self.objf_counting else: if objf_numpy: objf = self.objf_counting else: objf = self.objf_python_return_numpy if hessian_method.startswith('numdifftools'): import numdifftools as nd if hessian_method == 'python': h = hessian(objf, x, scalar=self.minimizing, perturbation=1e-4, zero_offset=self.jacobian_zero_offset, args=args) elif hessian_method == 'numdifftools_forward': h = nd.Hessian(objf, method='forward')(x) elif hessian_method == 'numdifftools_reverse': h = nd.Hessian(objf, method='reverse')(x) elif hessian_method == 'numdifftools_central': h = nd.Hessian(objf, method='central')(x) # Up the hessian fval count, set the fval back self.hess_fval_count += self.fval_iter - fval_iter self.fval_iter = fval_iter if return_numpy: return np.array(h) if type(h) is list else h else: return h if type(h) is list else h.tolist() def jacobian(self, x, *args): ''' jacobi - doesn't support jacobian_perturbation, jacobian_zero_offset, jacobian_order python - doesn't support jacobian_order ''' self.jac_iter += 1 fval_iter, jacobian_method, objf_numpy, jacobian_numpy = self.fval_iter, self.jacobian_method, self.objf_numpy, self.jacobian_numpy return_numpy = type(x) is not list if jacobian_method == 'analytical': if objf_numpy and not return_numpy: x = np.array(x) elif return_numpy and not objf_numpy: x = x.tolist() j = self.solver_analytical_jac(x, *args) else: # if the jacobian method doesn't speak numpy, convert x to a list if not jacobian_numpy: x = x if type(x) is list else x.tolist() else: x = np.array(x) # If the objf doesn't speak numpy but the jacobian does, use the converter if self.minimizing: if objf_numpy: objf = self.objf_numpy_minimizing else: objf = self.objf_python_minimizing else: if not jacobian_numpy: if objf_numpy: objf = self.objf_numpy_return_python else: objf = self.objf_counting else: if objf_numpy: objf = self.objf_counting else: objf = self.objf_python_return_numpy if jacobian_method.startswith('numdifftools'): import numdifftools as nd numdifftools_func = nd.Gradient if self.minimizing else nd.Jacobian step = self.jacobian_perturbation*x step[np.where(step==0)] = self.jacobian_zero_offset elif jacobian_method.startswith('jacobi'): from jacobi import jacobi if jacobian_method == 'python': j = jacobian(objf, x, scalar=self.minimizing, perturbation=self.jacobian_perturbation, zero_offset=self.jacobian_zero_offset, args=args) elif jacobian_method == 'numdifftools_forward': j = numdifftools_func(objf, method='forward', order=self.jacobian_order, step=step)(x) elif jacobian_method == 'numdifftools_reverse': j = numdifftools_func(objf, method='reverse', order=self.jacobian_order, step=step)(x) elif jacobian_method == 'numdifftools_central': j = numdifftools_func(objf, method='central', order=self.jacobian_order, step=step)(x) elif jacobian_method == 'jacobi_forward': j = jacobi(objf, x, method=1)[0] elif jacobian_method == 'jacobi_central': j = jacobi(objf, x, method=0)[0] elif jacobian_method == 'jacobi_backward': j = jacobi(objf, x, method=1)[0] # Up the jacobian fval count, set the fval back self.jac_fval_count += self.fval_iter - fval_iter self.fval_iter = fval_iter # Handle the return value - doesn't matter what type the method returns # if (not (jacobian_numpy ^ return_numpy)) or (return_numpy and jacobian_numpy): # return j # elif return_numpy and not jacobian_numpy: # return np.array(j) # elif jacobian_numpy and not return_numpy: # return j.tolist() if return_numpy: return np.array(j) if type(j) is list else j else: return j if type(j) is list else j.tolist() def solve(self, x0, args=()): self.fval_iter = self.jac_iter = self.jac_fval_count = self.hess_iter = self.hess_fval_count = 0 return_numpy = type(x0) is not list if self.solver_numpy: x0 = np.array(x0) elif return_numpy: x0 = x0.tolist() process_root = False method = self.method if method == 'newton_system': sln, niter = newton_system(self.objf, x0, jac=self.jacobian, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, damping=self.damping, solve_func=self.matrix_solver) elif method == 'newton_system_line_search': sln, niter = newton_system(self.objf, x0, jac=self.jacobian, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, damping=self.damping, line_search=True, solve_func=self.matrix_solver, require_progress=False, check_numbers=True, Armijo=True) elif method == 'newton_system_line_search_progress': sln, niter = newton_system(self.objf, x0, jac=self.jacobian, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, damping=self.damping, line_search=True, solve_func=self.matrix_solver, require_progress=True, check_numbers=True, Armijo=True) elif method == 'homotopy_solver': sln, niter = homotopy_solver(self.objf, x0, jac=self.jacobian, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, damping=self.damping, line_search=True, solve_func=self.matrix_solver) elif method == 'broyden2_python': sln, niter = broyden2(x0, self.objf, self.jacobian, xtol=self.xtol, maxiter=self.maxiter, args=args) elif method == 'fixed_point': sln, niter = fixed_point(residual_to_fixed_point(self.objf), x0, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, damping=self.damping) elif method == 'fixed_point_aitken': sln, niter = fixed_point_aitken(residual_to_fixed_point(self.objf), x0, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, damping=self.damping) elif method == 'fixed_point_gdem': sln, niter = fixed_point_gdem(residual_to_fixed_point(self.objf), x0, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, damping=self.damping) elif method == 'fixed_point_anderson': sln, niter = fixed_point_anderson_residual(residual_to_fixed_point(self.objf), x0, xtol=self.xtol, args=args, ytol=self.ytol, maxiter=self.maxiter, window_size=5, reg=1e-8, initial_iterations=5, acc_damping=0.3, damping=1.0) elif method in scipy_root_options_set: process_root = True jacobian_method = self.jacobian_method jac = self.jacobian if jacobian_method != 'scipy' else None result = root(self.objf, x0, args=args, method=method, jac=jac, tol=self.xtol) elif method in scipy_minimize_options_set: process_root = True jacobian_method = self.jacobian_method jac = self.jacobian if jacobian_method != 'scipy' else None hess = self.hessian if method in scipy_requires_hessian_options_set else None result = minimize(self.objf, x0, args=args, method=method, jac=jac, tol=self.xtol, hess=hess) if process_root: sln = result.x if not return_numpy: sln = sln.tolist() elif return_numpy: sln = np.array(sln) return sln 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, } # 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 # Make some dedicated references to the solvers which accet generic types e.g. mpmath # These should be used instead of e.g. newton anywhere mpmath/sympy/etc types get thrown in # This avoid numba issues newton_generic = newton secant_generic = secant bisect_generic = bisect brenth_generic = brenth