'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. Copyright (C) 2018 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, sublicense, 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 math import cos, erf, exp, isnan, log, pi, sin, sqrt import pytest from fluids.numerics import ( SolverInterface, array_as_tridiagonals, assert_close, assert_close1d, assert_close2d, best_bounding_bounds, chebder, chebint, chebval, chebval_ln_tau, chebval_ln_tau_and_der, chebval_ln_tau_and_der2, chebval_ln_tau_and_der3, cumsum, derivative, hessian, assert_close3d, exp_cheb, exp_cheb_and_der, exp_cheb_and_der2, exp_cheb_and_der3, exp_cheb_ln_tau, exp_cheb_ln_tau_and_der, exp_cheb_ln_tau_and_der2, fit_integral_linear_extrapolation, fit_integral_over_T_linear_extrapolation, full, horner, is_monotonic, is_poly_positive, isclose, jacobian, linspace, 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, min_max_ratios, newton_system, poly_fit_integral_over_T_value, poly_fit_integral_value, polyint, polyint_over_x, polylog2, polynomial_offset_scale, secant, sincos, fixed_point, solve_2_direct, solve_3_direct, solve_4_direct, solve_tridiagonal, std, subset_matrix, translate_bound_f_jac, translate_bound_func, translate_bound_jac, tridiagonals_as_array, trunc_exp_numpy, trunc_log_numpy, zeros, is_increasing, argsort1d, fixed_point_to_residual, residual_to_fixed_point, broyden2, fixed_point_aitken, fixed_point_gdem, fixed_point_anderson, ) from fluids.numerics import numpy as np assert_allclose = np.testing.assert_allclose def test_error_functions(): data = [1.0, 2.0, 3.0] calc = [.99, 2.01, 3.2] assert_close(max_abs_error(data, calc), 0.2, rtol=1e-13) assert_close(max_abs_rel_error(data, calc), 0.06666666666666672, rtol=1e-13) assert_close(max_squared_error(data, calc), 0.04000000000000007, rtol=1e-13) assert_close(max_squared_rel_error(data, calc), 0.004444444444444451, rtol=1e-13) assert_close(mean_abs_error(data, calc), 0.07333333333333332, rtol=1e-13) assert_close(mean_abs_rel_error(data, calc), 0.027222222222222207, rtol=1e-13) assert_close(mean_squared_error(data, calc), 0.013400000000000023, rtol=1e-13) assert_close(mean_squared_rel_error(data, calc), 0.0015231481481481502, rtol=1e-13) def test_sincos(): N = 10**1 for v in linspace(0.0, 2.0*pi, N): a, b = sincos(v) assert_close(a, sin(v), rtol=1e-14) assert_close(b, cos(v), rtol=1e-14) for v in linspace(-100.0, 100.0, N): a, b = sincos(v) assert_close(a, sin(v), rtol=1e-14) assert_close(b, cos(v), rtol=1e-14) def test_bisect_log_exp_terminations(): from math import exp, log from fluids.numerics import bisect def to_solve(x): try: return exp(x) except: return -1 assert 709.782712893384 == bisect(to_solve, 600, 800, xtol=1e-16) def to_solve(x): x = 10**x try: return log(x) except: return 1.0 assert -323.60724533877976 == bisect(to_solve, -300, -400, xtol=1e-16) def test_cumsum(): assert_close1d(cumsum([1,2,3,4,5]), [1, 3, 6, 10, 15]) assert_close1d(cumsum([1]), [1]) def test_interp(): from fluids.numerics import interp # Real world test data a = [0.29916, 0.29947, 0.31239, 0.31901, 0.32658, 0.33729, 0.34202, 0.34706, 0.35903, 0.36596, 0.37258, 0.38487, 0.38581, 0.40125, 0.40535, 0.41574, 0.42425, 0.43401, 0.44788, 0.45259, 0.47181, 0.47309, 0.49354, 0.49924, 0.51653, 0.5238, 0.53763, 0.54806, 0.55684, 0.57389, 0.58235, 0.59782, 0.60156, 0.62265, 0.62649, 0.64948, 0.65099, 0.6687, 0.67587, 0.68855, 0.69318, 0.70618, 0.71333, 0.72351, 0.74954, 0.74965] b = [0.164534, 0.164504, 0.163591, 0.163508, 0.163439, 0.162652, 0.162224, 0.161866, 0.161238, 0.160786, 0.160295, 0.15928, 0.159193, 0.157776, 0.157467, 0.156517, 0.155323, 0.153835, 0.151862, 0.151154, 0.14784, 0.147613, 0.144052, 0.14305, 0.140107, 0.138981, 0.136794, 0.134737, 0.132847, 0.129303, 0.127637, 0.124758, 0.124006, 0.119269, 0.118449, 0.113605, 0.113269, 0.108995, 0.107109, 0.103688, 0.102529, 0.099567, 0.097791, 0.095055, 0.087681, 0.087648] xs = np.linspace(0.29, 0.76, 100) ys = [interp(xi, a, b) for xi in xs.tolist()] ys_numpy = np.interp(xs, a, b) assert_allclose(ys, ys_numpy, atol=1e-12, rtol=1e-11) # Test custom extrapolation method xs = [1,2,3] ys = [.1, .2, .3] assert_close(interp(3.5, xs, ys, extrapolate=True), .35, rtol=1e-15) assert_close(interp(0, xs, ys, extrapolate=True), 0, rtol=1e-15) assert_close(interp(-1, xs, ys, extrapolate=True), -.1, rtol=1e-15) assert_close(interp(-100, xs, ys, extrapolate=True), -10, rtol=1e-15) assert_close(interp(10, xs, ys, extrapolate=True), 1, rtol=1e-15) assert_close(interp(10.0**30, xs, ys, extrapolate=True), 10.0**29, rtol=1e-15) def test_zeros(): # Only support up to 4d so far assert_allclose(zeros(2), np.zeros(2), atol=0) assert_allclose(zeros((2, )), np.zeros((2, )), atol=0) assert_allclose(zeros((2, 3)), np.zeros((2, 3)), atol=0) assert_allclose(zeros((2, 3, 4)), np.zeros((2, 3, 4)), atol=0) assert_allclose(zeros((2, 3, 4, 5)), np.zeros((2, 3, 4, 5)), atol=0) def test_full(): assert_allclose(full(2, 1), np.full(2, 1), atol=0) assert_allclose(full((2, ), 1), np.full((2, ), 1), atol=0) assert_allclose(full((2, 3), 1), np.full((2, 3), 1), atol=0) assert_allclose(full((2, 3, 4), 1), np.full((2, 3, 4), 1), atol=0) assert_allclose(full((2, 3, 4, 5), 1), np.full((2, 3, 4, 5), 1), atol=0) def test_splev(): from scipy.interpolate import splev from fluids.numerics import py_splev # Originally Dukler_XA_tck tck = [np.array([-2.4791105294648372, -2.4791105294648372, -2.4791105294648372, -2.4791105294648372, 0.14360803483759585, 1.7199938263676038, 1.7199938263676038, 1.7199938263676038, 1.7199938263676038]), np.array([0.21299880246561081, 0.16299733301915248, -0.042340970712679615, -1.9967836909384598, -2.9917366639619414, 0.0, 0.0, 0.0, 0.0]), 3] my_tck = [tck[0].tolist(), tck[1].tolist(), tck[2]] xs = np.linspace(-3, 2, 100).tolist() # test extrapolation ys_scipy = splev(xs, tck, ext=0) ys = [py_splev(xi, my_tck, ext=0) for xi in xs] assert_allclose(ys, ys_scipy) # test truncating to side values ys_scipy = splev(xs, tck, ext=3) ys = [py_splev(xi, my_tck, ext=3) for xi in xs] assert_allclose(ys, ys_scipy) # Test returning zeros for bad values ys_scipy = splev(xs, tck, ext=1) ys = [py_splev(xi, my_tck, ext=1) for xi in xs] assert_allclose(ys, ys_scipy) # Test raising an error when extrapolating is not allowed with pytest.raises(ValueError): py_splev(xs[0], my_tck, ext=2) with pytest.raises(ValueError): splev(xs[0], my_tck, ext=2) def test_bisplev(): from scipy.interpolate import bisplev from fluids.numerics import py_bisplev as my_bisplev tck = [np.array([0.0, 0.0, 0.0, 0.0, 0.0213694, 0.0552542, 0.144818, 0.347109, 0.743614, 0.743614, 0.743614, 0.743614]), np.array([0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0]), np.array([1.0001228445490002, 0.9988161050974387, 0.9987070557919563, 0.9979385859402731, 0.9970983069823832, 0.96602540121758, 0.955136014969614, 0.9476842472211648, 0.9351143114374392, 0.9059649602818451, 0.9218915266550902, 0.9086000082864022, 0.8934758292610783, 0.8737960765592091, 0.83185251064324, 0.8664296734965998, 0.8349705397843921, 0.809133298969704, 0.7752206120745123, 0.7344035693011536, 0.817047920445813, 0.7694560150930563, 0.7250979336267909, 0.6766754605968431, 0.629304180420512, 0.7137237030611423, 0.6408238328161417, 0.5772000233279148, 0.504889627280836, 0.440579886434288, 0.6239736474980684, 0.5273646894226224, 0.43995388722059986, 0.34359277007615313, 0.26986439252143746, 0.5640689738382749, 0.4540959882735219, 0.35278120580740957, 0.24364672351604122, 0.1606942128340308]), 3, 1] my_tck = [tck[0].tolist(), tck[1].tolist(), tck[2].tolist(), tck[3], tck[4]] xs = np.linspace(0, 1, 10) zs = np.linspace(0, 1, 10) ys_scipy = bisplev(xs, zs, tck) ys = my_bisplev(xs, zs, my_tck) assert_allclose(ys, ys_scipy) ys_scipy = bisplev(0.5, .7, tck) ys = my_bisplev(.5, .7, my_tck) assert_allclose(ys, ys_scipy) def test_linspace(): from fluids.numerics import linspace calc = linspace(-3,10, endpoint=True, num=8) expect = np.linspace(-3,10, endpoint=True, num=8) assert_allclose(calc, expect) calc = linspace(-3,10, endpoint=False, num=20) expect = np.linspace(-3,10, endpoint=False, num=20) assert_allclose(calc, expect) calc = linspace(0,1e-10, endpoint=False, num=3) expect = np.linspace(0,1e-10, endpoint=False, num=3) assert_allclose(calc, expect) calc = linspace(0,1e-10, endpoint=False, num=2) expect = np.linspace(0,1e-10, endpoint=False, num=2) assert_allclose(calc, expect) calc = linspace(0,1e-10, endpoint=False, num=1) expect = np.linspace(0,1e-10, endpoint=False, num=1) assert_allclose(calc, expect) calc, calc_step = linspace(0,1e-10, endpoint=False, num=2, retstep=True) expect, expect_step = np.linspace(0,1e-10, endpoint=False, num=2, retstep=True) assert_allclose(calc, expect) assert_allclose(calc_step, expect_step) calc, calc_step = linspace(0,1e-10, endpoint=False, num=1, retstep=True) expect, expect_step = np.linspace(0,1e-10, endpoint=False, num=1, retstep=True) assert_allclose(calc, expect) assert isnan(calc_step) # Cannot compare against numpy expect_step - it did not use to give nan in older versions calc, calc_step = linspace(100, 1000, endpoint=False, num=21, retstep=True) expect, expect_step = np.linspace(100, 1000, endpoint=False, num=21, retstep=True) assert_allclose(calc, expect) assert_allclose(calc_step, expect_step) def test_logspace(): from fluids.numerics import logspace calc = logspace(3,10, endpoint=True, num=8) expect = np.logspace(3,10, endpoint=True, num=8) assert_allclose(calc, expect) calc = logspace(3,10, endpoint=False, num=20) expect = np.logspace(3,10, endpoint=False, num=20) assert_allclose(calc, expect) calc = logspace(0,1e-10, endpoint=False, num=3) expect = np.logspace(0,1e-10, endpoint=False, num=3) assert_allclose(calc, expect) calc = logspace(0,1e-10, endpoint=False, num=2) expect = np.logspace(0,1e-10, endpoint=False, num=2) assert_allclose(calc, expect) calc = logspace(0,1e-10, endpoint=False, num=1) expect = np.logspace(0,1e-10, endpoint=False, num=1) assert_allclose(calc, expect) calc = logspace(0,1e-10, endpoint=False, num=2) expect = np.logspace(0,1e-10, endpoint=False, num=2) assert_allclose(calc, expect) calc = logspace(0,1e-10, endpoint=False, num=1) expect = np.logspace(0,1e-10, endpoint=False, num=1) assert_allclose(calc, expect) calc = logspace(100, 200, endpoint=False, num=21) expect = np.logspace(100, 200, endpoint=False, num=21) assert_allclose(calc, expect) def test_diff(): from fluids.numerics import diff test_arrs = [np.ones(10), np.zeros(10), np.arange(1, 10), np.arange(1, 10)*25.1241251, (np.arange(1, 10)**1.2), (10.1 + np.arange(1, 10)**20), (10.1 + np.linspace(-100, -10, 9)), (np.logspace(-10, -100, 19)**1.241), (np.logspace(10, 100, 15)**1.241) ] for test_arr in test_arrs: arr = test_arr.tolist() for n in range(5): diff_np = np.diff(arr, n=n) diff_py = diff(arr, n=n) assert_allclose(diff_np, diff_py) assert tuple(diff([1,2,3], n=0)) == (1,2,3) with pytest.raises(Exception): diff([1,2,3], n=-1) def test_fit_integral_linear_extrapolation(): coeffs = [-6.496329615255804e-23,2.1505678500404716e-19, -2.2204849352453665e-16, 1.7454757436517406e-14, 9.796496485269412e-11, -4.7671178529502835e-08, 8.384926355629239e-06, -0.0005955479316119903, 29.114778709934264] Tmin, Tmax = 50.0, 1000.0 Tmin_value, Tmax_value = 29.10061916353635, 32.697696220612684 Tmin_slope, Tmax_slope = 9.512557609301246e-06, 0.005910807286115391 int_coeffs = polyint(coeffs) T_int_T_coeffs, log_coeff = polyint_over_x(coeffs) def func(T): if T < Tmin: Cp = (T - Tmin)*Tmin_slope + Tmin_value elif T > Tmax: Cp = (T - Tmax)*Tmax_slope + Tmax_value else: Cp = horner(coeffs, T) return Cp assert_allclose(func(300), 29.12046448327871, rtol=1e-12) Ts = [0, 1, 25, 49, 50, 51, 500, 999, 1000, 1001, 2000, 50000] T_ends = [0, Tmin, Tmin*2.0, Tmax, Tmax*2.0] expect_values = [0.0, 0.0, 29.10014829193469, -29.10014829193469, 727.50656106565, -727.50656106565, 1425.9184530725483, -1425.9184530725483, 1455.0190674798057, -1455.0190674798057, 1484.119654811151, -1484.119654811151, 14588.023573849947, -14588.023573849947, 30106.09421115182, -30106.09421115182, 30138.789058157658, -30138.789058157658, 30171.489709781912, -30171.489709781912, 65791.88892182804, -65791.88892182804, 8728250.050849704, -8728250.050849704, -1455.0190674798057, 1455.0190674798057, -1425.918919187871, 1425.918919187871, -727.5125064141557, 727.5125064141557, -29.100614407257353, 29.100614407257353, 0.0, 0.0, 29.100587331345196, -29.100587331345196, 13133.004506370142, -13133.004506370142, 28651.075143672013, -28651.075143672013, 28683.76999067785, -28683.76999067785, 28716.470642302105, -28716.470642302105, 64336.869854348224, -64336.869854348224, 8726795.031782223, -8726795.031782223, -2910.0427925248396, 2910.0427925248396, -2880.942644232905, 2880.942644232905, -2182.53623145919, 2182.53623145919, -1484.1243394522915, 1484.1243394522915, -1455.023725045034, 1455.023725045034, -1425.923137713689, 1425.923137713689, 11677.980781325108, -11677.980781325108, 27196.05141862698, -27196.05141862698, 27228.746265632813, -27228.746265632813, 27261.446917257068, -27261.446917257068, 62881.84612930319, -62881.84612930319, 8725340.008057179, -8725340.008057179, -30138.789058157658, 30138.789058157658, -30109.688909865723, 30109.688909865723, -29411.282497092005, 29411.282497092005, -28712.870605085107, 28712.870605085107, -28683.76999067785, 28683.76999067785, -28654.669403346503, 28654.669403346503, -15550.765484307707, 15550.765484307707, -32.6948470058378, 32.6948470058378, 0.0, 0.0, 32.70065162425453, -32.70065162425453, 35653.09986367037, -35653.09986367037, 8698111.261791546, -8698111.261791546, -65791.88892182804, 65791.88892182804, -65762.7887735361, 65762.7887735361, -65064.38236076238, 65064.38236076238, -64365.97046875548, 64365.97046875548, -64336.869854348224, 64336.869854348224, -64307.769267016876, 64307.769267016876, -51203.86534797808, 51203.86534797808, -35685.794710676215, 35685.794710676215, -35653.09986367037, 35653.09986367037, -35620.39921204612, 35620.39921204612, 0.0, 0.0, 8662458.161927875, -8662458.161927875] numericals = [] analyticals = [] analyticals2 = [] for Tend in T_ends: for Tdiff in Ts: for (T1, T2) in zip([Tend, Tdiff], [Tdiff, Tend]): analytical = fit_integral_linear_extrapolation(T1, T2, int_coeffs, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope) analytical2 = (poly_fit_integral_value(T2, int_coeffs, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope) - poly_fit_integral_value(T1, int_coeffs, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope)) # numerical = quad(func, T1, T2, epsabs=1.49e-14, epsrel=1.49e-14)[0] # numericals.append(numerical) analyticals.append(analytical) analyticals2.append(analytical2) # print(analyticals) # assert_allclose(analyticals, numericals, rtol=1e-7) # assert_allclose(analyticals2, numericals, rtol=1e-7) assert_allclose(analyticals, expect_values, rtol=1e-12) assert_allclose(analyticals2, expect_values, rtol=1e-12) # Cannot have temperatures of 0 absolute for integrals over T cases Ts = [1e-9, 1, 25, 49, 50, 51, 500, 999, 1000, 1001, 2000, 50000] T_ends = [1e-9, Tmin, Tmin*2.0, Tmax, Tmax*2.0] expect_values = [0.0, 0.0, 603.0500198952521, -603.0500198952521, 696.719996723752, -696.719996723752, 716.3030057880156, -716.3030057880156, 716.890916983322, -716.890916983322, 717.467185070396, -717.467185070396, 783.9851859794122, -783.9851859794122, 805.3820356163964, -805.3820356163964, 805.4147468212567, -805.4147468212567, 805.4474311329551, -805.4474311329551, 829.8928106482913, -829.8928106482913, 1199.8352295965128, -1199.8352295965128, -716.890916983322, 716.890916983322, -113.84089708806982, 113.84089708806982, -20.170920259569826, 20.170920259569826, -0.5879111953062761, 0.5879111953062761, 0.0, 0.0, 0.5762680870740127, -0.5762680870740127, 67.09426899609028, -67.09426899609028, 88.4911186330744, -88.4911186330744, 88.5238298379347, -88.5238298379347, 88.5565141496331, -88.5565141496331, 113.00189366496929, -113.00189366496929, 482.94431261319096, -482.94431261319096, -737.0618021522312, 737.0618021522312, -134.01178225697902, 134.01178225697902, -40.34180542847902, 40.34180542847902, -20.75879636421547, 20.75879636421547, -20.170885168909194, 20.170885168909194, -19.59461708183518, 19.59461708183518, 46.923383827181084, -46.923383827181084, 68.3202334641652, -68.3202334641652, 68.3529446690255, -68.3529446690255, 68.38562898072391, -68.38562898072391, 92.8310084960601, -92.8310084960601, 462.77342744428177, -462.77342744428177, -805.4147468212567, 805.4147468212567, -202.36472692600452, 202.36472692600452, -108.69475009750452, 108.69475009750452, -89.11174103324097, 89.11174103324097, -88.5238298379347, 88.5238298379347, -87.94756175086069, 87.94756175086069, -21.42956084184442, 21.42956084184442, -0.03271120486030554, 0.03271120486030554, 0.0, 0.0, 0.03268431169840369, -0.03268431169840369, 24.47806382703459, -24.47806382703459, 394.42048277525623, -394.42048277525623, -829.8928106482913, 829.8928106482913, -226.8427907530391, 226.8427907530391, -133.17281392453913, 133.17281392453913, -113.58980486027556, 113.58980486027556, -113.00189366496929, 113.00189366496929, -112.42562557789527, 112.42562557789527, -45.90762466887901, 45.90762466887901, -24.510775031894894, 24.510775031894894, -24.47806382703459, 24.47806382703459, -24.445379515336185, 24.445379515336185, 0.0, 0.0, 369.94241894822164, -369.94241894822164] numericals = [] analyticals = [] analyticals2 = [] for Tend in T_ends: for Tdiff in Ts: for (T1, T2) in zip([Tend, Tdiff], [Tdiff, Tend]): analytical = fit_integral_over_T_linear_extrapolation(T1, T2, T_int_T_coeffs, log_coeff, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope) analytical2 = (poly_fit_integral_over_T_value(T2, T_int_T_coeffs, log_coeff, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope) - poly_fit_integral_over_T_value(T1, T_int_T_coeffs, log_coeff, Tmin, Tmax, Tmin_value, Tmax_value, Tmin_slope, Tmax_slope)) # numerical = quad(lambda T: func(float(T))/T, T1, T2, epsabs=1.49e-12, epsrel=1.49e-14)[0] # numericals.append(numerical) analyticals.append(analytical) analyticals2.append(analytical2) # # assert_allclose(analyticals, numericals, rtol=1e-7) # assert_allclose(analyticals2, numericals, rtol=1e-7) assert_allclose(analyticals, expect_values, rtol=1e-11) assert_allclose(analyticals2, expect_values, rtol=1e-11) def test_best_bounding_bounds(): def to_solve(x): return exp(x) - 100 vals = best_bounding_bounds(0, 5, to_solve, xs_pos=[4.831, 4.6054], ys_pos= [25.38, 0.0288], xs_neg=[4, 4.533, 4.6051690], ys_neg=[-45.40, -6.933, -0.0001139]) assert_allclose(vals, (4.605169, 4.6054, -0.0001139, 0.0288), rtol=1e-12) vals = best_bounding_bounds(4.60517018598, 5, to_solve, xs_pos=[4.831, 4.6054], ys_pos= [25.38, 0.0288], xs_neg=[4, 4.533, 4.6051690], ys_neg=[-45.40, -6.933, -0.0001139]) assert_allclose(vals, (4.60517018598, 4.6054, -8.091802783383173e-10, 0.0288), rtol=1e-12) vals = best_bounding_bounds(0, 4.60517018599, to_solve, xs_pos=[4.831, 4.6054], ys_pos= [25.38, 0.0288], xs_neg=[4, 4.533, 4.6051690], ys_neg=[-45.40, -6.933, -0.0001139]) assert_allclose(vals, (4.605169, 4.60517018599, -0.0001139, 1.908233571157325e-10), rtol=1e-12) vals = best_bounding_bounds(0, 4.60517018599, fa=to_solve(0), fb=to_solve(4.60517018599), xs_pos=[4.831, 4.6054], ys_pos= [25.38, 0.0288], xs_neg=[4, 4.533, 4.6051690], ys_neg=[-45.40, -6.933, -0.0001139]) assert_allclose(vals, (4.605169, 4.60517018599, -0.0001139, 1.908233571157325e-10), rtol=1e-12) def test_is_poly_positive(): assert not is_poly_positive([4, 3, 2, 1]) for high in range(0, 30, 5): assert is_poly_positive([4, 3, 2, 1], domain=(0, 10**high)) coeffs_4alpha = [2.1570803657937594e-10, 2.008831101045556e-06, -0.004656598178209313, 2.8575882247542514] assert not is_poly_positive(coeffs_4alpha) assert is_poly_positive(coeffs_4alpha, domain=(0, 511)) assert is_poly_positive(coeffs_4alpha, domain=(-10000, 511)) assert not is_poly_positive(coeffs_4alpha, domain=(-20000, 511)) assert not is_poly_positive(coeffs_4alpha, domain=(-15000, 511)) assert not is_poly_positive(coeffs_4alpha, domain=(-13000, 511)) assert not is_poly_positive(coeffs_4alpha, domain=(-11500, 511)) def test_translate_bound_func(): def rosen_test(x): x, y = x return (1.0 - x)**2 + 100.0*(y - x**2)**2 f, into, outof = translate_bound_func(rosen_test, low=[-2, -.2], high=[3.0, 4]) point = [.6, .7] in_exp = [0.0800427076735365, -1.2992829841302609] assert_allclose(into(point), in_exp, rtol=1e-12) assert_allclose(outof(in_exp), point, rtol=1e-12) assert f(into([1, 1])) < 1e-20 assert_allclose(outof(into([1, 1])), [1, 1], rtol=1e-12) f, into, outof = translate_bound_func(rosen_test, bounds=[[-2, 3], [-.2, 4]]) point = [.6, .7] in_exp = [0.0800427076735365, -1.2992829841302609] assert_allclose(into(point), in_exp, rtol=1e-12) assert_allclose(outof(in_exp), point, rtol=1e-12) assert f(into([1, 1])) < 1e-20 assert_allclose(outof(into([1, 1])), [1, 1], rtol=1e-12) def test_translate_bound_jac(): from scipy.optimize import rosen_der def rosen_test(x): x, y = x return (1.0 - x)**2 + 100.0*(y - x**2)**2 j, into, outof = translate_bound_jac(rosen_der, low=[-2, -.2], high=[3.0, 4]) f, into, outof = translate_bound_func(rosen_test, low=[-2, -.2], high=[3.0, 4]) point = [3, -2] jac_num = jacobian(f, point, perturbation=1e-8) jac_anal = j(point) assert_allclose(jac_num, jac_anal, rtol=1e-6) def test_translate_bound_f_jac(): from scipy.optimize import rosen_der def rosen_test(x): x, y = x return (1.0 - x)**2 + 100.0*(y - x**2)**2 low, high = [-2, -.2], [3.0, 4] f_j, into, outof = translate_bound_f_jac(rosen_test, rosen_der, low=low, high=high) point = [3, -2] f0, j0 = f_j(point) f0_check = translate_bound_func(rosen_test, low=low, high=high)[0](point) assert_allclose(f0_check, f0, rtol=1e-13) j0_check = translate_bound_jac(rosen_der, low=low, high=high)[0](point) assert_allclose(j0_check, j0, rtol=1e-13) def test_translate_bound_f_jac_multivariable(): def cons_f(x): return [x[0]**2 + x[1], x[0]**2 - x[1]] def cons_J(x): return[[2*x[0], 1], [2*x[0], -1]] new_f_j, translate_into, translate_outof = translate_bound_f_jac(cons_f, cons_J, bounds=[(-10, 10), (-5, 5)]) assert_close2d(new_f_j([-2.23, 3.45])[1], jacobian(lambda g: new_f_j(g, )[0], [-2.23, 3.45], scalar=False), rtol=3e-6) def test_solve_direct(): A = [[1.0,2.53252], [34.34, .5342]] B = [1.1241, .54354] assert_close1d(np.linalg.solve(A, B), solve_2_direct(A, B), rtol=1e-14) assert type(solve_2_direct(A, B)) is list A = [[1.0,2.53252, 54.54], [34.34, .5342, .545], [12.43, .545, .55555]] B = [1.1241, .54354, 1.22333] assert_close1d(np.linalg.solve(A, B), solve_3_direct(A, B), rtol=5e-14) assert type(solve_3_direct(A, B)) is list A = [[1.0,2.53252, 54.54, .235], [34.34, .5342, .545, .223], [12.43, .545, .55555, 33.33], [1.11, 2.2, 3.33, 4.44]] B = [1.1241, .54354, 1.22333, 9.009] ans = solve_4_direct(A, B) assert_close1d(np.linalg.solve(A, B), ans, rtol=5e-14) assert type(ans) is list def test_polylog2(): x = polylog2(0.5) assert_close(x, 0.5822405264516294) xs = linspace(0,0.99999, 50) # ys_act = [float(polylog(2, x)) for x in xs] # from sympy import polylog ys_act = [0.0, 0.020513035768572635, 0.0412401364927588, 0.06218738124039796, 0.08336114665629184, 0.10476812813791354, 0.12641536301777412, 0.1483102559926201, 0.1704606070746889, 0.19287464238138674, 0.21556104812821067, 0.23852900824703252, 0.261788246119884, 0.2853490709994786, 0.309222429784819, 0.33341996493707843, 0.3579540794622072, 0.38283801005840257, 0.4080859097363812, 0.43371294147827794, 0.45973538481992837, 0.4861707576383267, 0.5130379559237574, 0.540357414944675, 0.5681512960135646, 0.596443704089335, 0.6252609427828415, 0.6546318150738004, 0.6845879803506546, 0.7151643814663312, 0.7463997596771124, 0.7783372810645774, 0.811025306032588, 0.8445183447984893, 0.878878258148156, 0.9141757868110273, 0.9504925291123206, 0.9879235426949309, 1.026580835488432, 1.0665981582977615, 1.108137763432647, 1.1514002456979586, 1.1966394380910048, 1.244186068718536, 1.2944877067946645, 1.3481819162579485, 1.4062463083287482, 1.4703641942000052, 1.5441353484206717, 1.644808936992927] ys = [polylog2(x) for x in xs] assert_close1d(ys, ys_act, rtol=1E-7, atol=1E-10) def test_std(): inputs = ([1.0,5.0,11.0,4.0], [1.0,-5.0,11.0,-4.0], [1e12,-5.e13,11e14,-4e13], [1, 2, 3, 4], [-1, -2, -3, 4], [14, 8, 11, 10, 7, 9, 10, 11, 10, 15, 5, 10] ) for thing in inputs: assert_close(std(thing), np.std(thing), rtol=1e-14) def test_min_max_ratios(): actual = [1,2,3,4,5] calculated = [.9, 2.1, 3.05, 3.8, 5.5] min_ratio_np, max_ratio_np = np.min(np.array(calculated)/actual), np.max(np.array(calculated)/actual) assert_close1d([min_ratio_np, max_ratio_np], min_max_ratios(actual, calculated), rtol=1e-14) # Case with a zero match actual = [1,2,3,0,5] calculated = [.9, 2.1, 3.05, 0.0, 5.5] assert_close1d(min_max_ratios(actual, calculated), (0.9, 1.1), rtol=0) # Case with a zero mismatch actual = [1,2,3,0,5] calculated = [.9, 2.1, 3.05, 1, 5.5] assert_close1d(min_max_ratios(actual, calculated), (0.9, 10.0), rtol=0) def test_exp_cheb_fit_ln_tau(): coeffs = [-5.922664830406188, -3.6003367212635444, -0.0989717205896406, 0.05343895281736921, -0.02476759166597864, 0.010447569392539213, -0.004240542036664352, 0.0017273355647560718, -0.0007199858491173661, 0.00030714447101984343, -0.00013315510546685339, 5.832551964424226e-05, -2.5742454514671165e-05, 1.143577875153956e-05, -5.110008470393668e-06, 2.295229193177706e-06, -1.0355920205401548e-06, 4.690917226601865e-07, -2.1322112805921556e-07, 9.721709759435981e-08, -4.4448656630335925e-08, 2.0373327115630335e-08, -9.359475430792408e-09, 4.308620855930645e-09, -1.9872392620357004e-09, 9.181429297400179e-10, -4.2489342599871804e-10, 1.969051449668413e-10, -9.139573819982871e-11, 4.2452263926406886e-11, -1.9768853221080462e-11, 9.190537220149508e-12, -4.2949394041258415e-12, 1.9981863386142606e-12, -9.396025624219817e-13, 4.335282133283158e-13, -2.0410756418343112e-13, 1.0455525334407412e-13, -4.748978987834107e-14, 2.7630675525358583e-14] Tmin, Tmax, Tc = 233.22, 646.15, 647.096 expect = 0.031264474019763455 expect_d, expect_d2 = -0.00023379922039411865, -1.0453010755999069e-07 xmin, xmax = log(1-Tmin/Tc), log(1-Tmax/Tc) offset, scale = polynomial_offset_scale(xmin, xmax) coeffs_d = chebder(coeffs, m=1, scl=scale) coeffs_d2 = chebder(coeffs_d, m=1, scl=scale) coeffs_d3 = chebder(coeffs_d2, m=1, scl=scale) T = 500 calc = exp_cheb_ln_tau(T, Tc, coeffs, offset, scale) assert 0 == exp_cheb_ln_tau(700, Tc, coeffs, offset, scale) assert_close(expect, calc) calc2 = exp_cheb_ln_tau_and_der(T, Tc, coeffs, coeffs_d, offset, scale) assert (0,0) == exp_cheb_ln_tau_and_der(700, Tc, coeffs, coeffs_d, offset, scale) assert_close(expect, calc2[0]) assert_close(expect_d, calc2[1]) calc3 = exp_cheb_ln_tau_and_der2(T, Tc, coeffs, coeffs_d, coeffs_d2, offset, scale) assert (0,0,0) == exp_cheb_ln_tau_and_der2(700, Tc, coeffs, coeffs_d, coeffs_d2, offset, scale) assert_close(expect, calc3[0]) assert_close(expect_d, calc3[1]) assert_close(expect_d2, calc3[2]) def test_chebval_ln_tau(): Tmin, Tmax = 178.18, 591.0 Tc = 591.75 coeffs = [18231.740838720892, -18598.514785409734, 5237.841944302821, -1010.5549489362293, 147.88312821848922, -17.412144225239444, 1.7141064359038864, -0.14493639179363527, 0.01073811633477817, -0.0007078634084791702, 4.202655964036239e-05, -2.274648068123497e-06, 1.1239490049774759e-07] T = 500 xmin, xmax = log(1-Tmin/Tc), log(1-Tmax/Tc) offset, scale = polynomial_offset_scale(xmin, xmax) coeffs_d = chebder(coeffs, m=1, scl=scale) coeffs_d2 = chebder(coeffs_d, m=1, scl=scale) coeffs_d3 = chebder(coeffs_d2, m=1, scl=scale) calc = chebval_ln_tau(T, Tc, coeffs, offset, scale) assert 0 == chebval_ln_tau(600, Tc, coeffs, offset, scale) expect = 24498.131947622023 expect_d, expect_d2, expect_d3 = -100.77476795241955, -0.6838185834436981, -0.012093191904152178 assert_close(expect, calc) calc2 = chebval_ln_tau_and_der(T, Tc, coeffs, coeffs_d, offset, scale) assert (0,0) == chebval_ln_tau_and_der(600, Tc, coeffs, coeffs_d, offset, scale) assert_close(expect, calc2[0]) assert_close(expect_d, calc2[1]) calc3 = chebval_ln_tau_and_der2(T, Tc, coeffs, coeffs_d, coeffs_d2, offset, scale) assert (0,0,0) == chebval_ln_tau_and_der2(600, Tc, coeffs, coeffs_d, coeffs_d2, offset, scale) assert_close(expect, calc3[0]) assert_close(expect_d, calc3[1]) assert_close(expect_d2, calc3[2]) calc4 = chebval_ln_tau_and_der3(T, Tc, coeffs, coeffs_d, coeffs_d2, coeffs_d3, offset, scale) assert (0,0,0,0) == chebval_ln_tau_and_der3(600, Tc, coeffs, coeffs_d, coeffs_d2, coeffs_d3, offset, scale) assert_close(expect, calc4[0]) assert_close(expect_d, calc4[1]) assert_close(expect_d2, calc4[2]) assert_close(expect_d3, calc4[3]) def test_exp_cheb(): xmin, xmax = (309.0, 591.72) coeffs = [12.570668791524573, 3.1092695610681673, -0.5485217707981505, 0.11115875762247596, -0.01809803938553478, 0.003674911307077089, -0.00037626163070525465, 0.0001962813915017403, 6.120764548889213e-05, 3.602752453735203e-05] x = 400.0 offset, scale = polynomial_offset_scale(xmin, xmax) expect = 157186.81766860923 calc = exp_cheb(x, coeffs, offset, scale) assert_close(calc, expect, rtol=1e-14) coeffs_d = chebder(coeffs, m=1, scl=scale) coeffs_d2 = chebder(coeffs_d, m=1, scl=scale) coeffs_d3 = chebder(coeffs_d2, m=1, scl=scale) der_num = derivative(exp_cheb, x, args=(coeffs, offset, scale), dx=x*1e-7) der_analytical = exp_cheb_and_der(x, coeffs, coeffs_d, offset, scale)[1] assert_close(der_num, der_analytical, rtol=1e-7) assert_close(der_analytical, 4056.277312107932, rtol=1e-14) der_num = derivative(lambda *args: exp_cheb_and_der(*args)[1], x, args=(coeffs, coeffs_d, offset, scale), dx=x*1e-7) der_analytical = exp_cheb_and_der2(x, coeffs, coeffs_d, coeffs_d2, offset, scale)[-1] assert_close(der_analytical, 81.34302144188977, rtol=1e-14) assert_close(der_num, der_analytical, rtol=1e-7) der_num = derivative(lambda *args: exp_cheb_and_der2(*args)[-1], x, args=(coeffs, coeffs_d, coeffs_d2, offset, scale), dx=x*1e-7) der_analytical = exp_cheb_and_der3(x, coeffs, coeffs_d, coeffs_d2, coeffs_d3, offset, scale)[-1] assert_close(der_num, der_analytical, rtol=1e-7) assert_close(der_analytical, 1.105438780935656, rtol=1e-14) vals = exp_cheb_and_der3(x, coeffs, coeffs_d, coeffs_d2, coeffs_d3, offset, scale) assert_close1d(vals, (157186.81766860923, 4056.277312107932, 81.34302144188977, 1.105438780935656), rtol=1e-14) vals = exp_cheb_and_der2(x, coeffs, coeffs_d, coeffs_d2, offset, scale) assert_close1d(vals, (157186.81766860923, 4056.277312107932, 81.34302144188977), rtol=1e-14) vals = exp_cheb_and_der(x, coeffs, coeffs_d, offset, scale) assert_close1d(vals, (157186.81766860923, 4056.277312107932), rtol=1e-14) def test_cheb(): Tmin, Tmax = 50, 1500.0 toluene_TRC_cheb_fit = [194.9993931442641, 135.143566535142, -31.391834328585, -0.03951841213554952, 5.633110876073714, -3.686554783541794, 1.3108038668007862, -0.09053861376310801, -0.2614279887767278, 0.24832452742026911, -0.15919652548841812, 0.09374295717647019, -0.06233192560577938, 0.050814520356653126, -0.046331125185531064, 0.0424579816955023, -0.03739513702085129, 0.031402017733109244, -0.025212485578021915, 0.01939423141593144, -0.014231480849538403, 0.009801281575488097, -0.006075456686871594, 0.0029909809015365996, -0.0004841890018462136, -0.0014991199985455728, 0.0030051480117581075, -0.004076901418829215, 0.004758297389532928, -0.005096275567543218, 0.00514099984344718, -0.004944736724873944, 0.004560044671604424, -0.004037777783658769, 0.0034252408915679267, -0.002764690626354871, 0.0020922734527478726, -0.0014374230267101273, 0.0008226963858916081, -0.00026400260413972365, -0.0002288377348015347, 0.0006512726893767029, -0.0010030137199867895, 0.0012869214641443305, -0.001507857723972772, 0.001671575150882565, -0.0017837100581746812, 0.001848935469520696, -0.0009351605848800237] toluene_TRC_cheb_fit_copy = [v for v in toluene_TRC_cheb_fit] offset, scale = polynomial_offset_scale(Tmin, Tmax) val = chebval(300, toluene_TRC_cheb_fit, offset, scale) assert_close(val, 104.46956642594124, rtol=1e-14) d1_coeffs = chebder(toluene_TRC_cheb_fit, m=1, scl=scale) d2_coeffs = chebder(toluene_TRC_cheb_fit, m=2, scl=scale) d2_2_coeffs = chebder(d1_coeffs, m=1, scl=scale) val = chebval(300, d1_coeffs, offset, scale) assert_close(val, 0.36241217517888635, rtol=1e-14) val = chebval(300, d2_coeffs, offset, scale) assert_close(val, -6.445511348110282e-06, rtol=1e-14) val = chebval(300, d2_2_coeffs, offset, scale) assert_close(val, -6.445511348110282e-06, rtol=1e-14) assert d2_2_coeffs == d2_coeffs int_coeffs = chebint(toluene_TRC_cheb_fit, m=1, lbnd=0, scl=1/scale) assert_close(chebval(300, int_coeffs, offset, scale), -83708.18079449862, rtol=1e-10) assert toluene_TRC_cheb_fit == toluene_TRC_cheb_fit_copy def test_cheb_more(): c = [1, 2, 3] x = 0.5 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 0.5 assert_close(result, expected, rtol=1e-13) c = [1] x = 1.0 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 1.0 assert_close(result, expected, rtol=1e-13) c = [1, 2] x = -1.0 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = -1.0 assert_close(result, expected, rtol=1e-13) c = [0, 0, 1] x = 2.0 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 7.0 assert_close(result, expected, rtol=1e-13) c = [1, -1, 1, -1, 1] x = 0.1 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 1.1368 assert_close(result, expected, rtol=1e-13) c = [0.1, 0.2, 0.3, 0.4, 0.5] x = 3.0 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 333.8999999999999 assert_close(result, expected, rtol=1e-13) c = [10, 20, 30, 40, 50] x = -5.0 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 222030.0 assert_close(result, expected, rtol=1e-13) c = [1e-05, 2e-05, 3e-05, 4e-05] x = 100.0 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 160.58998000000003 assert_close(result, expected, rtol=1e-13) c = [1, 0, 0, 0, 1] x = 3.141592653589793 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 702.3158930633044 assert_close(result, expected, rtol=1e-13) c = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] x = 0.75 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 6.9560546875 assert_close(result, expected, rtol=1e-13) c = [1, -1, 1, -1, 1, -1, 1, -1] x = -0.5 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 1.5 assert_close(result, expected, rtol=1e-13) c = [1.0, 0.5, 0.3333333333333333, 0.25, 0.2, 0.16666666666666666] x = 1.5 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 20.116666666666664 assert_close(result, expected, rtol=1e-13) c = [2.718281828459045, 3.141592653589793, 2.718, 3.142] x = 0.25 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = -1.0346950081435065 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = -1 scale = 2 result = chebval(x, c, offset, scale) expected = -2.0 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] x = 1.0 offset = 10 scale = 100 result = chebval(x, c, offset, scale) expected = 5877283503.0 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 0.001 scale = 1 result = chebval(x, c, offset, scale) expected = 0.508006 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 1e-06 scale = 1 result = chebval(x, c, offset, scale) expected = 0.5000080000060003 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 1000.0 scale = 1 result = chebval(x, c, offset, scale) expected = 6008000.5 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 1000000.0 scale = 1 result = chebval(x, c, offset, scale) expected = 6000008000000.5 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 1000000000.0 scale = 1 result = chebval(x, c, offset, scale) expected = 6.000000008e+18 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 0 scale = 0.001 result = chebval(x, c, offset, scale) expected = -1.9989985 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 0 scale = 1e-06 result = chebval(x, c, offset, scale) expected = -1.9999989999985 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 0 scale = 1000.0 result = chebval(x, c, offset, scale) expected = 1500998.0 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 0 scale = 1000000.0 result = chebval(x, c, offset, scale) expected = 1500000999998.0 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 0 scale = 1000000000.0 result = chebval(x, c, offset, scale) expected = 1.500000001e+18 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 1e-06 scale = 1e-06 result = chebval(x, c, offset, scale) expected = -1.9999969999865 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 1000000.0 scale = 1000000.0 result = chebval(x, c, offset, scale) expected = 13500002999998.0 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = -1000000000.0 scale = 2000000000.0 result = chebval(x, c, offset, scale) expected = -2.0 assert_close(result, expected, rtol=1e-13) c = [1e-09, 2e-09, 3e-09] x = 1000000000.0 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 6000000002.0 assert_close(result, expected, rtol=1e-13) c = [1000000000.0, 2000000000.0, 3000000000.0] x = 1e-09 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = -1999999998.0 assert_close(result, expected, rtol=1e-13) c = [42] x = 0.5 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 42.0 assert_close(result, expected, rtol=1e-13) c = [42] x = 1000 offset = 1000000.0 scale = 1e-06 result = chebval(x, c, offset, scale) expected = 42.0 assert_close(result, expected, rtol=1e-13) c = [42] x = -1000 offset = -1000000.0 scale = 1e-06 result = chebval(x, c, offset, scale) expected = 42.0 assert_close(result, expected, rtol=1e-13) c = [1e-09] x = 1000000000.0 offset = 1e-09 scale = 1000000000.0 result = chebval(x, c, offset, scale) expected = 1.0e-9 assert_close(result, expected, rtol=1e-13) c = [1000000000.0] x = 1e-09 offset = 1000000000.0 scale = 1e-09 result = chebval(x, c, offset, scale) expected = 1000000000.0 assert_close(result, expected, rtol=1e-13) c = [1, 2] x = 0.5 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 2.0 assert_close(result, expected, rtol=1e-13) c = [1, 2] x = 1000 offset = 1000000.0 scale = 1e-06 result = chebval(x, c, offset, scale) expected = 2000001.002 assert_close(result, expected, rtol=1e-13) c = [1, 2] x = -1000 offset = -1000000.0 scale = 1e-06 result = chebval(x, c, offset, scale) expected = -1999999.002 assert_close(result, expected, rtol=1e-13) c = [1e-09, 2e-09] x = 1000000000.0 offset = 1e-09 scale = 1000000000.0 result = chebval(x, c, offset, scale) expected = 2000000000.0 assert_close(result, expected, rtol=1e-13) c = [1000000000.0, 2000000000.0] x = 1e-09 offset = 1000000000.0 scale = 1e-09 result = chebval(x, c, offset, scale) expected = 2.000000001e+18 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 0.5 offset = 0 scale = 1 result = chebval(x, c, offset, scale) expected = 0.5 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = 1000 offset = 1000000.0 scale = 1e-06 result = chebval(x, c, offset, scale) expected = 6000002011998.0 assert_close(result, expected, rtol=1e-13) c = [1, 2, 3] x = -1000 offset = -1000000.0 scale = 1e-06 result = chebval(x, c, offset, scale) expected = 5999998011998.0 assert_close(result, expected, rtol=1e-13) c = [1e-09, 2e-09, 3e-09] x = 1000000000.0 offset = 1e-09 scale = 1000000000.0 result = chebval(x, c, offset, scale) expected = 6.0e+27 assert_close(result, expected, rtol=1e-13) c = [1000000000.0, 2000000000.0, 3000000000.0] x = 1e-09 offset = 1000000000.0 scale = 1e-09 result = chebval(x, c, offset, scale) expected = 6.000000002e+27 assert_close(result, expected, rtol=1e-13) c = [0] x = 3.141592653589793 offset = -1000000.0 scale = 1000000.0 result = chebval(x, c, offset, scale) expected = 0.0 assert_close(result, expected, rtol=1e-13) c = [1, 0] x = 1e-09 offset = 1000000000.0 scale = 1e-09 result = chebval(x, c, offset, scale) expected = 1.0 assert_close(result, expected, rtol=1e-13) c = [0, 0, 1] x = 1000000000.0 offset = -1000000000.0 scale = 2000000000.0 result = chebval(x, c, offset, scale) expected = 7.999999992e+36 assert_close(result, expected, rtol=1e-13) c = [1e-15] x = 1000000000000000.0 offset = 1000000000000000.0 scale = 1e-15 result = chebval(x, c, offset, scale) expected = 1.0e-15 assert_close(result, expected, rtol=1e-13) c = [1000000000000000.0, 1e-15] x = 1e-15 offset = 1e-15 scale = 1000000000000000.0 result = chebval(x, c, offset, scale) expected = 1.0e+15 assert_close(result, expected, rtol=1e-13) def test_chebder_more(): c = [1, 2, 3, 4, 5] m = 0 scl = 1.0 result = chebder(c, m, scl) expected = [1, 2, 3, 4, 5] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 1 scl = 1.0 result = chebder(c, m, scl) expected = [14.0, 52.0, 24, 40] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 1 scl = 2.0 result = chebder(c, m, scl) expected = [28.0, 104.0, 48.0, 80.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 1 scl = 0.5 result = chebder(c, m, scl) expected = [7.0, 26.0, 12.0, 20.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 2 scl = 1.0 result = chebder(c, m, scl) expected = [172.0, 96.0, 240] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 2 scl = 2.0 result = chebder(c, m, scl) expected = [688.0, 384.0, 960.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 2 scl = 0.5 result = chebder(c, m, scl) expected = [43.0, 24.0, 60.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 3 scl = 1.0 result = chebder(c, m, scl) expected = [96.0, 960.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 3 scl = 2.0 result = chebder(c, m, scl) expected = [768.0, 7680.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 3 scl = 0.5 result = chebder(c, m, scl) expected = [12.0, 120.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 4 scl = 1.0 result = chebder(c, m, scl) expected = [960.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 4 scl = 2.0 result = chebder(c, m, scl) expected = [15360.0] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3, 4, 5] m = 4 scl = 0.5 result = chebder(c, m, scl) expected = [60.0] assert_close1d(result, expected, rtol=1e-13) c = [1] m = 1 scl = 1.0 result = chebder(c, m, scl) expected = [] assert_close1d(result, expected, rtol=1e-13) c = [1, 2] m = 1 scl = 1.0 result = chebder(c, m, scl) expected = [2] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3] m = 3 scl = 1.0 result = chebder(c, m, scl) expected = [] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3] m = 4 scl = 1.0 result = chebder(c, m, scl) expected = [] assert_close1d(result, expected, rtol=1e-13) c = [1, 2, 3] m = 5 scl = 1.0 result = chebder(c, m, scl) expected = [] assert_close1d(result, expected, rtol=1e-13) c = [1, 0, 1, 0, 1] m = 2 scl = 1.0 result = chebder(c, m, scl) expected = [36.0, 0.0, 48] assert_close1d(result, expected, rtol=1e-13) c = [1e-05, 2e-05, 3e-05, 4e-05] m = 1 scl = 100000.0 result = chebder(c, m, scl) expected = [14.0, 12.0, 24.0] assert_close1d(result, expected, rtol=1e-13) c = [100000.0, 200000.0, 300000.0, 400000.0] m = 1 scl = 1e-05 result = chebder(c, m, scl) expected = [14.0, 12.000000000000002, 24.0] assert_close1d(result, expected, rtol=1e-13) def test_chebint_more(): c = [1, 2, 3, 4, 5] m = 0 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [1.0, 2.0, 3.0, 4.0, 5.0] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 1 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [-1.0, -0.5, -0.5, -0.33333333333333337, 0.5, 0.5] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 1 lbnd = 1 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.3333333333333335, -0.5, -0.5, -0.33333333333333337, 0.5, 0.5] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 1 lbnd = -1 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [-0.33333333333333337, -0.5, -0.5, -0.33333333333333337, 0.5, 0.5] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 1 lbnd = 0 scl = 2.0 result = chebint(c, m, lbnd, scl) expected = [-2.0, -1.0, -1.0, -0.6666666666666667, 1.0, 1.0] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 1 lbnd = 0 scl = 0.5 result = chebint(c, m, lbnd, scl) expected = [-0.5, -0.25, -0.25, -0.16666666666666669, 0.25, 0.25] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 2 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.10416666666666669, -0.75, -0.04166666666666666, -0.16666666666666666, -0.10416666666666667, 0.05, 0.041666666666666664] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 2 lbnd = 1 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [-0.36250000000000016, 0.5833333333333335, -0.04166666666666666, -0.16666666666666666, -0.10416666666666667, 0.05, 0.041666666666666664] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 2 lbnd = -1 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [-0.09583333333333335, -0.08333333333333337, -0.04166666666666666, -0.16666666666666666, -0.10416666666666667, 0.05, 0.041666666666666664] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 2 lbnd = 0 scl = 2.0 result = chebint(c, m, lbnd, scl) expected = [0.41666666666666674, -3.0, -0.16666666666666663, -0.6666666666666666, -0.4166666666666667, 0.2, 0.16666666666666666] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 3 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [-0.11458333333333334, 0.125, -0.14583333333333334, 0.010416666666666668, -0.027083333333333334, -0.014583333333333334, 0.004166666666666667, 0.002976190476190476] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 3 lbnd = 1 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.17827380952380967, -0.34166666666666684, 0.18750000000000003, 0.010416666666666668, -0.027083333333333334, -0.014583333333333334, 0.004166666666666667, 0.002976190476190476] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 3 lbnd = -1 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [-0.07410714285714287, -0.07500000000000002, 0.020833333333333322, 0.010416666666666668, -0.027083333333333334, -0.014583333333333334, 0.004166666666666667, 0.002976190476190476] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 4 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.023871527777777776, -0.04166666666666667, 0.028645833333333332, -0.019791666666666666, 0.003125, -0.003125, -0.0014632936507936508, 0.0002976190476190476, 0.00018601190476190475] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4, 5] m = 4 lbnd = 1 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [-0.03128720238095246, 0.08452380952380965, -0.08802083333333338, 0.035763888888888894, 0.003125, -0.003125, -0.0014632936507936508, 0.0002976190476190476, 0.00018601190476190475] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1] m = 1 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.0, 1.0] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2] m = 1 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.5, 1.0, 0.5] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [0] m = 1 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.0] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3] m = 5 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.0008680555555555565, -0.009114583333333336, 0.0013020833333333343, -0.004427083333333334, 0.0005208333333333334, -0.0007812499999999999, 8.680555555555556e-05, 3.720238095238095e-05] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 0, 1, 0, 1] m = 2 lbnd = 0 scl = 1.0 result = chebint(c, m, lbnd, scl) expected = [0.14583333333333334, 0.0, 0.125, 0.0, -0.0125, 0.0, 0.008333333333333333] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1e-05, 2e-05, 3e-05, 4e-05] m = 1 lbnd = 0 scl = 100000.0 result = chebint(c, m, lbnd, scl) expected = [-1.0, -0.5, -0.5, 0.5, 0.5] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [100000.0, 200000.0, 300000.0, 400000.0] m = 1 lbnd = 0 scl = 1e-05 result = chebint(c, m, lbnd, scl) expected = [-1.0, -0.5000000000000002, -0.5, 0.5000000000000001, 0.5] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3] m = 1 lbnd = 2 scl = 3 result = chebint(c, m, lbnd, scl) expected = [-46.5, -1.5, 1.5, 1.5] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) c = [1, 2, 3, 4] m = 2 lbnd = -1 scl = 0.5 result = chebint(c, m, lbnd, scl) expected = [0.08020833333333333, 0.0625, -0.0625, -0.041666666666666664, 0.015625, 0.0125] assert_close1d(result, expected, rtol=1e-13, atol=1e-13) def test_is_monotonic(): assert is_monotonic([1,2,3]) assert is_monotonic([3, 2, 1]) assert is_monotonic([1,1,2,3]) assert is_monotonic([3,3, 2, 1]) assert is_monotonic([1]) assert not is_monotonic([3,3,5, 2, 1]) assert not is_monotonic([1,2,3,1]) assert is_monotonic([-1, -2]) assert not is_monotonic([-1, -2, -3, -2]) assert is_monotonic([-2, -1, 0]) assert not is_monotonic([-2, -1, 0, -1]) def test_secant_cases_tough(): def to_solve(x): err = -0.1 if x < 0.5 else x - 0.6 # Bad case, completely flat, cannot convergewhen x < 0.5 but difficulties elsewhere too return err assert_close(secant(to_solve, 1e20, xtol=1e-12, ytol=1e-20, require_xtol=False, bisection=True, maxiter=200), 0.6) def to_solve(x): # Have to try multiple initial guesses to get a working point err = x**20.0 - 1.0 return err assert_close(secant(to_solve, .104, xtol=1e-12, ytol=1e-20, require_xtol=False, bisection=True, maxiter=1000, additional_guesses=True), 1) def to_solve(x): err = 10.0**14*(1.0*x**7 - 7.0*x**6 + 21.0*x**5 - 35.0*x**4 + 35.0*x**3 - 21.0*x**2 + 7.0*x - 1.0) return err assert_close(secant(to_solve, 1.013, xtol=1e-12, ytol=1e-20, require_xtol=False, bisection=True, maxiter=1000, additional_guesses=True), 0.9926199643387525) # Go through a set of cases collected online def newton_baffler(x): y = x - 6.25 if y < -0.25: return 0.75*y- 0.3125 elif y < 0.25: return 2.0*y else: return 0.75*y + 0.3125 def newton_baffler_d(x): # others are zero y = x - 6.25 if y < -0.25: return 0.75 elif y < 0.25: return 2.0 else: return 0.75 def flat_stanley(x): if x == 1: return 0 else: if x < 1: return -exp(log(1000) + log(1.0 - x) - 1.0/(x - 1.0)**2) return exp(log(1000) + log(x - 1.0) - 1.0/(x - 1.0)**2) # factor = (-1.0 if x < 1.0 else 1.0) # return factor*exp(log(1000) + log(abs(x - 1.0)) - 1.0/(x - 1.0)**2) def flat_stanley_d(x): if x == 1: return 0 else: if x < 1: return 1000*exp(-1.0/(x - 1.0)**2) + 2.0*(1000*x - 1000.0)*exp(-1.0/(x - 1.0)**2)/(x - 1.0)**3 return 1000*exp(-1.0/(x - 1.0)**2) + 2.0*(1000*x - 1000.0)*exp(-1.0/(x - 1.0)**2)/(x - 1.0)**3 def flat_stanley_d2(x): if x == 1: return 0 else: if x < 1: return (-(6.0 - 4.0/(x - 1.0)**2)*(1000*x - 1000.0)/(x - 1.0) + 4000.0)*exp(-1.0/(x - 1.0)**2)/(x - 1.0)**3 return (-(6.0 - 4.0/(x - 1.0)**2)*(1000*x - 1000.0)/(x - 1.0) + 4000.0)*exp(-1.0/(x - 1.0)**2)/(x - 1.0)**3 def flat_stanley_d3(x): if x == 1: return 0 else: if x < 1: return (-18000.0 + (1000*x - 1000.0)*(24.0 - 36.0/(x - 1.0)**2 + 8.0/(x - 1.0)**4)/(x - 1.0) + 12000.0/(x - 1.0)**2)*exp(-1.0/(x - 1.0)**2)/(x - 1.0)**4 return (-18000.0 + (1000*x - 1000.0)*(24.0 - 36.0/(x - 1.0)**2 + 8.0/(x - 1.0)**4)/(x - 1.0) + 12000.0/(x - 1.0)**2)*exp(-1.0/(x - 1.0)**2)/(x - 1.0)**4 def newton_pathological(x): if x == 0.0: return 0.0 else: factor = (-1.0 if x < 0.0 else 1.0) return factor*abs(x)**(1.0/3.0)*exp(-x**2) solver_1d_test_cases = [(lambda x: sin(x) - x/2, lambda x: cos(x) - 1/2, lambda x: -sin(x), lambda x: -cos(x), [0.1], None, None, 'TEST_ZERO #1'), (lambda x: 2*x - exp(-x), lambda x: 2 + exp(-x), lambda x: -exp(-x), lambda x: exp(-x), [1.], None, None, 'TEST_ZERO #2'), (lambda x: x*exp(-x), lambda x: -x*exp(-x) + exp(-x), lambda x: (x - 2)*exp(-x), lambda x: (3 - x)*exp(-x), [1.], None, None, 'TEST_ZERO #3'), (lambda x: exp(x) - 1.0/(10.0*x)**2, lambda x: exp(x) + 0.02/x**3, lambda x: exp(x) - 0.06/x**4, lambda x: exp(x) + 0.24/x**5, [1.], None, None, 'TEST_ZERO #4'), (lambda x: (x+3)*(x-1)**2, lambda x: (x - 1)**2 + (x + 3)*(2*x - 2), lambda x: 2*(3*x + 1), lambda x: 6, [1.], None, None, 'TEST_ZERO #5'), (lambda x: exp(x) - 2 - 1/(10*x)**2 + 2/(100*x)**3, lambda x: exp(x) + 1/(50*x**3) - 3/(500000*x**4), lambda x: exp(x) - 3/(50*x**4) + 3/(125000*x**5), lambda x: exp(x) + 6/(25*x**5) - 3/(25000*x**6), [1.], None, None, 'TEST_ZERO #6'), (lambda x: x**3, lambda x: 3*x**2, lambda x: 6*x, lambda x: 6.0, [1.], None, None, 'TEST_ZERO #7'), (lambda x: cos(x) - x, lambda x: -sin(x) - 1, lambda x: -cos(x), lambda x: sin(x), [1.], None, None, 'TEST_ZERO #8'), (newton_baffler, newton_baffler_d, lambda x: 0, lambda x: 0, [6.25 + 5.0, 6.25 - 1.0, 6.25 + 0.1], None, None, 'TEST_ZERO #9/Newton Baffler'), (lambda x: 20*x/(100*x**2 + 1), lambda x: -4000*x**2/(100*x**2 + 1)**2 + 20/(100*x**2 + 1), lambda x: 4000*x*(400*x**2/(100*x**2 + 1) - 3)/(100*x**2 + 1)**2, lambda x: 12000*(-400*x**2*(200*x**2/(100*x**2 + 1) - 1)/(100*x**2 + 1) + 400*x**2/(100*x**2 + 1) - 1)/(100*x**2 + 1)**2, [1., -0.14, 0.041], None, None, 'TEST_ZERO #10/Repeller'), (lambda x: (16 - x**4)/(16*x**4 + 0.00001), lambda x: -x**3*(16 - x**4)/(4*(x**4 + 6.25e-7)**2) - 4*x**3/(16*x**4 + 1.0e-5), lambda x: x**2*(2*x**4/(x**4 + 6.25e-7)**2 - (x**4 - 16)*(8*x**4/(x**4 + 6.25e-7) - 3)/(4*(x**4 + 6.25e-7)**2) - 12/(16*x**4 + 1.0e-5)), lambda x: 3*x*(-x**4*(8*x**4/(x**4 + 6.25e-7) - 3)/(x**4 + 6.25e-7)**2 + 3*x**4/(x**4 + 6.25e-7)**2 + (x**4 - 16)*(16*x**8/(x**4 + 6.25e-7)**2 - 12*x**4/(x**4 + 6.25e-7) + 1)/(2*(x**4 + 6.25e-7)**2) - 8/(16*x**4 + 1.0e-5)), [.25, 5., 1.1], None, None, 'TEST_ZERO #11/Pinhead'), (flat_stanley, flat_stanley_d, flat_stanley_d2, flat_stanley_d3, [2, 0.5, 4], None, None, 'TEST_ZERO #12/Flat Stanley'), (lambda x: 0.00000000001*(x - 100),lambda x: 1.00000000000000e-11, lambda x: 0, lambda x: 0, [100., 1.], None, None, 'TEST_ZERO #13/Lazy Boy'), (lambda x: 1.0/((x - 0.3)**2 + 0.01) + 1.0/((x - 0.9)**2 + 0.04) + 2.0*x - 5.2, lambda x: 1.0*(0.6 - 2*x)/((x - 0.3)**2 + 0.01)**2 + 1.0*(1.8 - 2*x)/((x - 0.9)**2 + 0.04)**2 + 2.0, lambda x: 1.0*(2*x - 1.8)*(4*x - 3.6)/((x - 0.9)**2 + 0.04)**3 + 1.0*(2*x - 0.6)*(4*x - 1.2)/((x - 0.3)**2 + 0.01)**3 - 2.0/((x - 0.3)**2 + 0.01)**2 - 2.0/((x - 0.9)**2 + 0.04)**2, lambda x: -1.0*(2*x - 1.8)*(4*x - 3.6)*(6*x - 5.4)/((x - 0.9)**2 + 0.04)**4 - 1.0*(2*x - 0.6)*(4*x - 1.2)*(6*x - 1.8)/((x - 0.3)**2 + 0.01)**4 + (8.0*x - 7.2)/((x - 0.9)**2 + 0.04)**3 + (8.0*x - 2.4)/((x - 0.3)**2 + 0.01)**3 + (16.0*x - 14.4)/((x - 0.9)**2 + 0.04)**3 + (16.0*x - 4.8)/((x - 0.3)**2 + 0.01)**3, [3., -0.5, 0, 2.12742], None, None, 'TEST_ZERO #14/Camel'), # (newton_pathological, [0.01, -0.25], None, None, 'TEST_ZERO #15/Newton Pathological'), (lambda x: pi*(x - 5.0)/180.0 - 0.8*sin(pi*x/180), lambda x: -0.00444444444444444*pi*cos(pi*x/180) + 0.00555555555555556*pi, lambda x: 2.46913580246914e-5*pi**2*sin(pi*x/180), lambda x: 1.37174211248285e-7*pi**3*cos(pi*x/180), [0., 5+180., 5.], None, None, 'TEST_ZERO #16/Kepler'), (lambda x: x**3 - 2*x - 5,lambda x: 3*x**2 - 2, lambda x: 6*x, lambda x: 6, [2., 3.], None, None, 'TEST_ZERO #17/Wallis example'), (lambda x: 10.0**14*(1.0*x**7 - 7.0*x**6 + 21.0*x**5 - 35.0*x**4 + 35.0*x**3 - 21.0*x**2 + 7.0*x - 1.0), lambda x: 700000000000000.0*x**6 - 4.2e+15*x**5 + 1.05e+16*x**4 - 1.4e+16*x**3 + 1.05e+16*x**2 - 4.2e+15*x + 700000000000000.0, lambda x: 4.2e+15*x**5 - 2.1e+16*x**4 + 4.2e+16*x**3 - 4.2e+16*x**2 + 2.1e+16*x - 4.2e+15, lambda x: 2.1e+16*x**4 - 8.4e+16*x**3 + 1.26e+17*x**2 - 8.4e+16*x + 2.1e+16, [.99, 1.013], None, None, 'TEST_ZERO #18'), (lambda x: cos(100.0*x) - 4.0*erf(30*x - 10), lambda x: -100.0*sin(100.0*x) - 240.0*exp(-(30*x - 10)**2)/sqrt(pi), lambda x: (432000.0*x - 144000.0)*exp(-100*(3*x - 1)**2)/sqrt(pi) - 10000.0*cos(100.0*x), lambda x: -86400000.0*(3*x - 1)**2*exp(-100*(3*x - 1)**2)/sqrt(pi) + 1000000.0*sin(100.0*x) + 432000.0*exp(-100*(3*x - 1)**2)/sqrt(pi), [0., 1., 0.5], None, None, 'TEST_ZERO #19'), # # From scipy (lambda x: x**2 - 2*x - 1, lambda x: 2*x - 2, lambda x: 2, lambda x: 0, [3,], None, None, 'Scipy/x**2 - 2*x - 1'), (lambda x: exp(x) - cos(x), lambda x: exp(x) + sin(x), lambda x: exp(x) + cos(x), lambda x: exp(x) - sin(x), [3,], None, None, 'Scipy/exp(x) - cos(x)'), (lambda x: x - 0.1, lambda x: 1, lambda x: 0, lambda x: 0, [-1e8, 1e7], None, None, 'Scipy/GH5555'), (lambda x: -0.1 if x < 0.5 else x - 0.6, lambda x: 0.0 if x < 0.5 else 1.0, lambda x: 0, lambda x: 0, [1e20, 1.], None, None, 'Scipy/GH5557'), # Fail at 0 is expected - 0 slope (lambda x: (x - 100.0)**2, lambda x: 2.0*x - 200.0, lambda x: 2., lambda x: 0, [10*(200.0 - 6.828499381469512e-06) / (2.0 + 6.828499381469512e-06)], None, None, 'Scipy/zero_der_nz_dp'), (lambda x: x**3 - x**2, lambda x: 3*x**2 - 2*x, lambda x: 2*(3*x - 1), lambda x: 6, [0., 0.5], None, None, 'Scipy/GH8904'), (lambda x: x**(1.00/9.0) - 9**(1.0/9), lambda x: 1/(9*x**(8/9)), lambda x: -8/(81*x**(17/9)), lambda x: 136/(729*x**(26/9)), [0.1], None, None, 'Scipy/GH8881'), # # From scipy optimization suite - root functions only (lambda x: sin(x) + sin(10.0 / 3.0 * x), lambda x: cos(x) + 10*cos(10*x/3)/3, lambda x: -(sin(x) + 100*sin(10*x/3)/9), lambda x: -(cos(x) + 1000*cos(10*x/3)/27), [2.7048, 3.18, 4.14, 4.62, 5.1, 5.58, 6.06, 7.02, 7.4952], 2.7, 7.5, 'Scipy/Problem02'), (lambda x: -sum(k * sin((k + 1) * x + k) for k in range(1, 6)), lambda x: -2*cos(2*x + 1) - 6*cos(3*x + 2) - 12*cos(4*x + 3) - 20*cos(5*x + 4) - 30*cos(6*x + 5), lambda x: 2*(2*sin(2*x + 1) + 9*sin(3*x + 2) + 24*sin(4*x + 3) + 50*sin(5*x + 4) + 90*sin(6*x + 5)), lambda x: 2*(4*cos(2*x + 1) + 27*cos(3*x + 2) + 96*cos(4*x + 3) + 250*cos(5*x + 4) + 540*cos(6*x + 5)), [-9.98, -8.0, -4.0, -2.0, 0.0, 2.0, 4, 8.0, 9.98], -10, 10, 'Scipy/Problem03'), (lambda x: -(1.4 - 3 * x) * sin(18.0 * x), lambda x: 18.0*(3*x - 1.4)*cos(18.0*x) + 3*sin(18.0*x), lambda x: -324.0*(3*x - 1.4)*sin(18.0*x) + 108.0*cos(18.0*x), lambda x: -(5832.0*(3*x - 1.4)*cos(18.0*x) + 2916.0*sin(18.0*x)), [0.0012, 0.12, 0.36, 0.48, 0.6, 0.72, 0.84, 1.08, 1.1988], 0, 1.2, 'Scipy/Problem05'), (lambda x: -(x + sin(x)) * exp(-x ** 2.0), lambda x: -2*x*(-x - sin(x))*exp(-x**2) + (-cos(x) - 1)*exp(-x**2), lambda x: (4*x*(cos(x) + 1) - 2*(x + sin(x))*(2*x**2 - 1) + sin(x))*exp(-x**2), lambda x: (4*x*(x + sin(x))*(2*x**2 - 3) - 6*x*sin(x) - 6*(2*x**2 - 1)*(cos(x) + 1) + cos(x))*exp(-x**2), [-9.98, -8.0, -4.0, -2.0, 0.0, 2.0, 4, 8.0, 9.98], -10, 10, 'Scipy/Problem06'), (lambda x: sin(x) + sin(10.0 / 3.0 * x) + log(x) - 0.84 * x + 3, lambda x: cos(x) + 10*cos(10*x/3)/3 - 21/25 + 1/x, lambda x: -(sin(x) + 100*sin(10*x/3)/9 + x**(-2)), lambda x: -cos(x) - 1000*cos(10*x/3)/27 + 2/x**3, [2.7048, 3.18, 4.14, 4.62, 5.1, 5.58, 6.06, 7.02, 7.4952], 2.7, 7.5, 'Scipy/Problem07'), (lambda x: -sum(k * cos((k + 1) * x + k) for k in range(1, 6)), lambda x: 2*sin(2*x + 1) + 6*sin(3*x + 2) + 12*sin(4*x + 3) + 20*sin(5*x + 4) + 30*sin(6*x + 5), lambda x: 2*(2*cos(2*x + 1) + 9*cos(3*x + 2) + 24*cos(4*x + 3) + 50*cos(5*x + 4) + 90*cos(6*x + 5)), lambda x: -2*(4*sin(2*x + 1) + 27*sin(3*x + 2) + 96*sin(4*x + 3) + 250*sin(5*x + 4) + 540*sin(6*x + 5)), [-9.98, -8.0, -4.0, -2.0, 0.0, 2.0, 4, 8.0, 9.98], -10, 10, 'Scipy/Problem08'), (lambda x: sin(x) + sin(2.0 / 3.0 * x), lambda x: 2*cos(2*x/3)/3 + cos(x), lambda x: -(4*sin(2*x/3)/9 + sin(x)), lambda x: -(8*cos(2*x/3)/27 + cos(x)), [3.1173, 4.83, 8.29, 10.02, 11.75, 13.48, 15.21, 18.67, 20.3827], 3.1, 20.4, 'Scipy/Problem09'), (lambda x: -x * sin(x), lambda x: -x*cos(x) - sin(x), lambda x: x*sin(x) - 2*cos(x), lambda x: x*cos(x) + 3*sin(x), [0.01, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 9.0, 9.99], 0, 10, 'Scipy/Problem10'), (lambda x: 2 * cos(x) + cos(2 * x), lambda x: -2*sin(x) - 2*sin(2*x), lambda x: -2*(cos(x) + 2*cos(2*x)), lambda x: 2*(sin(x) + 4*sin(2*x)), [-1.5629423451609221, -0.7853981633974483, 0.7853981633974483, 1.5707963267948966, 2.356194490192345, 3.141592653589793, 3.9269908169872414, 5.497787143782138, 6.275331325545612], -pi / 2, 2 * pi, 'Scipy/Problem11'), (lambda x: (sin(x))**3.0 + (cos(x))**3.0, lambda x: 3*sin(x)**2*cos(x) - 3*sin(x)*cos(x)**2, lambda x: 3*(-sin(x)**3 + 2*sin(x)**2*cos(x) + 2*sin(x)*cos(x)**2 - cos(x)**3), lambda x: 3*(-2*sin(x)**3 - 7*sin(x)**2*cos(x) + 7*sin(x)*cos(x)**2 + 2*cos(x)**3), [0.006283185307179587, 0.6283185307179586, 1.8849555921538759, 2.5132741228718345, 3.141592653589793, 3.7699111843077517, 4.39822971502571, 5.654866776461628, 6.276902121872407], 0, 2*pi, 'Scipy/Problem12'), (lambda x: -exp(-x) * sin(2.0 * pi * x), lambda x: exp(-x)*sin(2*pi*x) - 2*pi*exp(-x)*cos(2*pi*x), lambda x: (-sin(2*pi*x) + 4*pi**2*sin(2*pi*x) + 4*pi*cos(2*pi*x))*exp(-x), lambda x: (-12*pi**2*sin(2*pi*x) + sin(2*pi*x) - 6*pi*cos(2*pi*x) + 8*pi**3*cos(2*pi*x))*exp(-x), [0.004, 0.4, 1.2, 1.6, 2.0, 2.4, 2.8, 3.6, 3.996], 0, 4, 'Scipy/Problem14'), (lambda x: -(x - sin(x)) * exp(-x ** 2.0), lambda x: -2*x*(-x + sin(x))*exp(-x**2) + (cos(x) - 1)*exp(-x**2), lambda x: -(4*x*(cos(x) - 1) + 2*(x - sin(x))*(2*x**2 - 1) + sin(x))*exp(-x**2), lambda x: (4*x*(x - sin(x))*(2*x**2 - 3) + 6*x*sin(x) + 6*(2*x**2 - 1)*(cos(x) - 1) - cos(x))*exp(-x**2), [-9.98, -8.0, -4.0, -2.0, 0.0, 2.0, 4, 8.0, 9.98], -10, 10, 'Scipy/Problem20'), (lambda x: x * sin(x) + x * cos(2.0 * x), lambda x: -2.0*x*sin(2.0*x) + x*cos(x) + sin(x) + cos(2.0*x), lambda x: -x*sin(x) - 4.0*x*cos(2.0*x) - 4.0*sin(2.0*x) + 2*cos(x), lambda x: 8.0*x*sin(2.0*x) - x*cos(x) - 3*sin(x) - 12.0*cos(2.0*x), [0.01, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 9.0, 9.99], 0, 10, 'Scipy/Problem21'), (lambda x: exp(-3.0 * x) - (sin(x)) ** 3.0, lambda x: -3*sin(x)**2*cos(x) - 3*exp(-3*x), lambda x: 3*(sin(x)**3 - 2*sin(x)*cos(x)**2 + 3*exp(-3*x)), lambda x: 3*(7*sin(x)**2*cos(x) - 2*cos(x)**3 - 9*exp(-3*x)), [0.02, 2.0, 6.0, 8.0, 10.0, 12.0, 14.0, 18.0, 19.98], 0, 20, 'Scipy/Problem22'), # # from gsl (lambda x: sin(x), lambda x: cos(x), lambda x: -sin(x), lambda x: -cos(x), [3, 4, -4, -3, -1/3., 1], None, None, 'sin(x)'), (lambda x: cos(x), lambda x: -sin(x), lambda x: -cos(x), lambda x: sin(x), [0, 3, -3, 0], None, None, 'cos(x)'), (lambda x: x**20 - 1,lambda x: 20*x**19, lambda x: 380*x**18, lambda x: 6840*x**17, [0.1, 2], None, None, 'x^20 - 1'), # Numerical derivative at .1 is 0, failing secant # (lambda x: np.sign(x)*abs(x)**0.5, [-1.0/3, 1], None, None, 'sign(x)*sqrt(abs(x))'), (lambda x: x**2 - 1e-8, lambda x: 2*x, lambda x: 2, lambda x: 0, [0, 1], None, None, 'x^2 - 1e-8'), (lambda x: (x-1.0)**7, lambda x: 7*(x - 1.0)**6, lambda x: 42*(x - 1.0)**5, lambda x: 210*(x - 1.0)**4, [0.9995, 1.0002], None, None, '(x-1.0)**7'), # From Roots.jl # (lambda x: abs(x - 0.0), [0, 1], None, None, 'abs(x - 0.0)'), (lambda x: 1024*x**11 - 2816*x**9 + 2816*x**7 - 1232*x**5 + 220*x**3 - 11*x, lambda x: 11264*x**10 - 25344*x**8 + 19712*x**6 - 6160*x**4 + 660*x**2 - 11, lambda x: 88*x*(1280*x**8 - 2304*x**6 + 1344*x**4 - 280*x**2 + 15), lambda x: 264*(3840*x**8 - 5376*x**6 + 2240*x**4 - 280*x**2 + 5), [0, -1, 1, 21], None, None, '1024*x**11 - 2816*x**9 + 2816*x**7 - 1232*x**5 + 220*x**3 - 11*x'), (lambda x: 512*x**9 - 1024*x**7 + 672*x**5 - 160*x**3 +10*x, lambda x: 4608*x**8 - 7168*x**6 + 3360*x**4 - 480*x**2 + 10, lambda x: 192*x*(192*x**6 - 224*x**4 + 70*x**2 - 5), lambda x: 192*(1344*x**6 - 1120*x**4 + 210*x**2 - 5), [0, -1, 1, 7], None, None, '512*x**9 - 1024*x**7 + 672*x**5 - 160*x**3 +10*x') ] def test_secant_cases_internet(): solvers = [secant] solver_kwargs = [{'xtol': 1e-14, 'bisection': True, 'require_xtol': False, 'ytol': 1e-12, 'additional_guesses': True, 'require_eval': True}, ] names = ['fluids secant'] def print_err(fun): def f(x): err = fun(x) # print([err, x]) return err return f for solver, kwargs, name in zip(solvers, solver_kwargs, names): passes = 0 fails = 0 failed_fs = set() for f, fder, fder2, fder3, xs, low, high, note in solver_1d_test_cases: for x0 in xs: v = solver(print_err(f), x0, maxiter=3000, low=low, high=high, **kwargs) try: v = solver(print_err(f), x0, maxiter=3000, low=low, high=high, **kwargs) assert abs(f(v)) < 1e-10 passes += 1 except Exception as e: failed_fs.add(f) fails += 1 assert fails == 0 @pytest.mark.filterwarnings("ignore:divide by zero encountered in scalar divide") @pytest.mark.filterwarnings("ignore:invalid value encountered in power") @pytest.mark.filterwarnings("ignore:invalid value encountered in scalar power") def test_secant_cases_nan_inf(): import numpy as np from fluids.numerics import UnconvergedError def div_by_zero(x): err = 1 / (np.array(x) - 1) return err def sqrt_of_negative(x): x = np.array(x) return (x - 2) ** 0.5 with pytest.raises((ValueError, UnconvergedError)): secant(div_by_zero, 1) with pytest.raises((ValueError, UnconvergedError)): secant(div_by_zero, x0=2, x1=1) with pytest.raises((ValueError, UnconvergedError)): secant(sqrt_of_negative, 1) with pytest.raises((ValueError, UnconvergedError)): secant(sqrt_of_negative, x0=4, x1=1) # test we can solve :) assert secant(sqrt_of_negative, 2, x1=2.1) == 2 # Force the solver to go into the nan regime via a high damping factor with pytest.raises((ValueError, UnconvergedError)): secant(sqrt_of_negative, 3, x1=2.1, damping=4) def test_newton_halley_cases_internet(): debug = False fcalls = 0 from math import inf, isinf, isnan from fluids.numerics import newton def print_err(fun): def f(x): nonlocal fcalls fcalls += 1 try: err = fun(x) except: return inf if err.imag != 0.0: return inf # if debug: # print(f"x={x}, err={err} fder={fder(x)}") return err return f def print_fder(fun): def f(x): try: fder = fun(x) except: return inf if fder.imag != 0.0: return inf return fder return f solvers = [newton] solver_kwargs = [{'xtol': 1e-14, 'bisection': True, 'ytol': 1e-11, 'require_xtol': False, 'require_eval': True, 'additional_guesses': True}, ] names = ['fluids newton'] for solver, kwargs, name in zip(solvers, solver_kwargs, names): passes = 0 fails = 0 failed_fs = set() for f, fder, fder2, fder3, xs, low, high, note in solver_1d_test_cases: for x0 in xs: # if debug: # print(f'Solving with x0={x0}') try: # low=low, high=high v = solver(print_err(f), x0, print_fder(fder), low=low, high=high, maxiter=3000, **kwargs) v = solver(print_err(f), x0, print_fder(fder), low=low, high=high, fprime2=fder2, maxiter=3000, **kwargs) v = solver(print_err(f), x0, print_fder(fder), fprime2=fder2, maxiter=3000, **kwargs) assert not isinf(v) assert not isnan(v) assert abs(f(v)) < 1e-10 passes += 1 except Exception as e: # print('Failed', e, note, x0, low, high) failed_fs.add(f) fails += 1 assert fails == 0 def test_basic_newton_system(): # Define system of equations def system(inputs): x, y = inputs # Example system: # f1 = x^2 + y^2 - 4 = 0 # f2 = exp(x) - y = 0 return [x**2 + y**2 - 4, exp(x) - y] # Define Jacobian matrix def jacobian(inputs): x, y = inputs return [ [2*x, 2*y], [exp(x), -1] ] # Initial guess x0 = [1.0, 1.0] # Solve system solution, iterations = newton_system( f=system, x0=x0, jac=jacobian, xtol=1e-10, maxiter=100 ) # Check results assert iterations > 0 assert iterations < 10 # Should converge in about 6 iters # Verify solution satisfies equations final_residuals = system(solution) assert_close1d(final_residuals, [0, 0], atol=1e-6) def test_newton_system_no_iterations(): def test_objf_direct(inputs): x, T = inputs k = 0.12*exp(12581*(T-298.)/(298.*T)) return [120*x-75*k*(1-x), -x*(873-T)-11.0*(T-300)] def test_jac_not_called(inputs): raise ValueError("Should not be called") expect = [0.05995136780143791, 296.85996516970505] found, iterations = newton_system(test_objf_direct, expect, jac=test_jac_not_called, ytol=1e-7, xtol=None) assert iterations == 0 assert found == expect def test_newton_system_with_damping_function(): # can get up to all sorts of fun with damping, useful for very weird cases def damping_func_max_one_second(x, step, damping, *args): new_step = [] for v in step: if abs(v) > damping: if v < 0: new_step.append(-damping) else: new_step.append(damping) else: new_step.append(v) xnew = [xi+s for xi, s in zip(x, new_step)] return xnew Ts = [] def test_objf_with_damping(inputs): x, T = inputs Ts.append(T) k = 0.12*exp(12581*(T-298.)/(298.*T)) return [120*x-75*k*(1-x), -x*(873-T)-11.0*(T-300)] def test_jac_with_damping(inputs): x, T = inputs ans = [[9.0*exp(0.00335570469798658*(12581*T - 3749138.0)/T) + 120, (42.2181208053691/T - 0.00335570469798658*(12581*T - 3749138.0)/T**2)*(9.0*x - 9.0)*exp(0.00335570469798658*(12581*T - 3749138.0)/T),], [T - 873, x - 11.0]] return ans near_solution = [0.05995136780143791, 300] ans, iterations = newton_system(test_objf_with_damping, near_solution, jac=test_jac_with_damping, line_search=True, damping=1, ytol=1e-7, xtol=None, damping_func=damping_func_max_one_second) assert Ts[0:4] == [300, 299, 298, 297] Ts = [] ans, iterations = newton_system(test_objf_with_damping, near_solution, jac=test_jac_with_damping, line_search=True, damping=2, ytol=1e-7, xtol=None, damping_func=damping_func_max_one_second) assert [300, 298] == Ts[0:2] def test_newton_system_nan_inf_inputs_outputs(): for bad in (float("nan"), float("inf")): Ts = [] def test_objf_with_nan_iterations(inputs): x, T = inputs Ts.append(T) k = 0.12*exp(12581*(T-298.)/(298.*T)) if len(Ts) == 3: return [bad, bad] return [120*x-75*k*(1-x), -x*(873-T)-11.0*(T-300)] def test_jac_with_nan_point(inputs): x, T = inputs ans = [[9.0*exp(0.00335570469798658*(12581*T - 3749138.0)/T) + 120, (42.2181208053691/T - 0.00335570469798658*(12581*T - 3749138.0)/T**2)*(9.0*x - 9.0)*exp(0.00335570469798658*(12581*T - 3749138.0)/T),], [T - 873, x - 11.0]] return ans near_solution = [0.05995136780143791, 300] with pytest.raises(ValueError): newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point, line_search=False, check_numbers=True, ytol=1e-7, xtol=None) Ts = [] # The line search will allow the solver to try again ans, iterations = newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point, line_search=True, check_numbers=True, xtol=1e-12) assert_close1d(ans, [0.059951367801437914, 296.85996516970505]) # the lack of line search causes a failure for progress in (True, False): with pytest.raises(ValueError): Ts = [] ans, iterations = newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point, line_search=False, check_numbers=True, xtol=1e-12, require_progress=progress) def test_newton_jacobian_has_nan_inf(): Ts = [] def test_objf_with_nan_iterations(inputs): x, T = inputs Ts.append(T) k = 0.12*exp(12581*(T-298.)/(298.*T)) return [120*x-75*k*(1-x), -x*(873-T)-11.0*(T-300)] def test_jac_with_nan_point(inputs): x, T = inputs ans = [[9.0*exp(0.00335570469798658*(12581*T - 3749138.0)/T) + 120, (42.2181208053691/T - 0.00335570469798658*(12581*T - 3749138.0)/T**2)*(9.0*x - 9.0)*exp(0.00335570469798658*(12581*T - 3749138.0)/T),], [T - 873, x - 11.0]] # if len(Ts) == 3: ans[0][0] = float('inf') return ans near_solution = [0.05995136780143791, 300] with pytest.raises(ValueError): newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point, line_search=False, check_numbers=True, ytol=1e-7, xtol=None) Ts = [] # # The line search will allow the solver to try again ans, iterations = newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point, line_search=True, check_numbers=True, xtol=1e-12, jac_error_allowed=True) assert_close1d(ans, [0.059951367801437914, 296.85996516970505]) # Should also work without a linesearch Ts = [] ans, iterations = newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point, line_search=False, check_numbers=True, xtol=1e-12, jac_error_allowed=True) def test_jac_with_nan_point_first_only(inputs): x, T = inputs ans = [[9.0*exp(0.00335570469798658*(12581*T - 3749138.0)/T) + 120, (42.2181208053691/T - 0.00335570469798658*(12581*T - 3749138.0)/T**2)*(9.0*x - 9.0)*exp(0.00335570469798658*(12581*T - 3749138.0)/T),], [T - 873, x - 11.0]] if len(Ts) == 1: ans[0][0] = float('inf') return ans # Check that if the initial jacobian is bad, it gets caught Ts = [] with pytest.raises(ValueError): ans, iterations = newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point_first_only, line_search=True, check_numbers=True, xtol=1e-12, jac_error_allowed=False) Ts = [] # check that the initial jacobian won't stop the problem if we say yes ans, iterations = newton_system(test_objf_with_nan_iterations, near_solution, jac=test_jac_with_nan_point_first_only, line_search=False, check_numbers=True, xtol=1e-12, jac_error_allowed=True) assert_close1d(ans, [0.059951367801437914, 296.85996516970505]) def to_solve_newton_python(inputs): x, T = inputs if not isinstance(inputs, list): raise ValueError("Bad input type") k = 0.12*exp(12581*(T-298.)/(298.*T)) return [120*x-75*k*(1-x), -x*(873-T)-11.0*(T-300)] def to_solve_newton_numpy(inputs): if not isinstance(inputs, np.ndarray): raise ValueError("Bad input type") x, T = inputs k = 0.12*exp(12581*(T-298.)/(298.*T)) return np.array([120*x-75*k*(1-x), -x*(873-T)-11.0*(T-300)]) def to_solve_jac_newton_python(inputs): # derived with sympy if not isinstance(inputs, list): raise ValueError("Bad input type") x, T = inputs ans = [[9.0*exp(0.00335570469798658*(12581*T - 3749138.0)/T) + 120, (42.2181208053691/T - 0.00335570469798658*(12581*T - 3749138.0)/T**2)*(9.0*x - 9.0)*exp(0.00335570469798658*(12581*T - 3749138.0)/T),], [T - 873, x - 11.0]] return ans def to_solve_jac_newton_numpy(inputs): # derived with sympy if not isinstance(inputs, np.ndarray): raise ValueError("Bad input type") x, T = inputs ans = [[9.0*exp(0.00335570469798658*(12581*T - 3749138.0)/T) + 120, (42.2181208053691/T - 0.00335570469798658*(12581*T - 3749138.0)/T**2)*(9.0*x - 9.0)*exp(0.00335570469798658*(12581*T - 3749138.0)/T),], [T - 873, x - 11.0]] return np.array(ans) try: import warnings with warnings.catch_warnings(): warnings.simplefilter("ignore") try: import numdifftools as nd jacob_methods = ['analytical', 'python', 'numdifftools_forward'] except: jacob_methods = ['analytical', 'python'] try: import jacobi jacob_methods += ['jacobi_forward'] except: pass except: pass @pytest.mark.parametrize("jacob_method", jacob_methods) @pytest.mark.filterwarnings("ignore:Method") def test_SolverInterface_basics(jacob_method): solver = SolverInterface(method='newton_system', objf=to_solve_newton_python, jacobian_method=jacob_method, jac=to_solve_jac_newton_python) # Not testing convergence so start quite close to the point x0 = [.06, 297.] res = solver.solve(x0) expect = [0.05995136780143791, 296.85996516970505] assert_close1d(res, expect) res = solver.solve(np.array(x0)) assert isinstance(res, np.ndarray) j0 = jacobian(to_solve_newton_python, x0, scalar=False) assert_close2d(to_solve_jac_newton_python(x0), j0, rtol=1e-5) assert_close2d(to_solve_jac_newton_numpy(np.array(x0)), j0, rtol=1e-5) if jacob_method == 'python': working_methods = ('newton_system_line_search', 'homotopy_solver', 'hybr', 'lm', 'krylov', 'df-sane', 'Nelder-Mead', 'Powell', 'BFGS', 'Newton-CG', 'TNC', 'SLSQP') else: working_methods = ('newton_system_line_search', 'hybr', 'Powell', 'BFGS') for method in working_methods: solver = SolverInterface(method=method, objf=to_solve_newton_python, jacobian_method=jacob_method, jac=to_solve_jac_newton_python) res = solver.solve(x0) assert_close1d(res, expect, rtol=1e-5) assert type(res) is list if method in ('newton_system_line_search', 'hybr', 'Powell', 'BFGS'): res = solver.solve(np.array(x0)) assert isinstance(res, np.ndarray) j = solver.jacobian(np.array(x0)) assert isinstance(j, np.ndarray) for method in working_methods: solver = SolverInterface(method=method, objf=to_solve_newton_numpy, jacobian_method=jacob_method, objf_numpy=True, jac=to_solve_jac_newton_numpy) res = solver.solve(np.array(x0)) assert_close1d(res, expect, rtol=1e-5) assert isinstance(res, np.ndarray) if method in ('newton_system_line_search', 'hybr', 'Powell', 'BFGS'): res = solver.solve([1, 400.]) assert type(res) is list def fixed_point_1_func(x): conversion = 0.8 stoichiometry = [-1, 0.5, 2, 1, 0.1, 0.001, 0, 0] fixed_point_1_feed = [20, 10, 0, 0, 0, 0, 30, 15] recycle = [x_i for x_i in x] reactor_feed = [recycle_i + feed_i for recycle_i, feed_i in zip(recycle, fixed_point_1_feed)] effluent = [reactor_feed_i + (reactor_feed[0] * stoich_i * conversion) for reactor_feed_i, stoich_i in zip(reactor_feed, stoichiometry)] product = [effluent_i * 0.1 for effluent_i in effluent] thing = [effluent_i - product_i for effluent_i, product_i in zip(effluent, product)] return [xi- ti for xi, ti in zip(x, thing)] fixed_point_1_guess = [20, 10, 0, 0, 0, 0, 30, 15] fixed_point_1_expect = [4.390243902439022, 177.80487804878027, 351.2195121951214, 175.6097560975607, 17.560975609756067, 0.17560975609756074, 269.9999999999998, 134.9999999999999] def test_fixed_point_process(): result, iterations = fixed_point(fixed_point_1_func, fixed_point_1_guess, xtol=1e-9, maxiter=1000) assert iterations < 250 # 235 last time assert_close1d(result, fixed_point_1_expect, rtol=1e-9) # Test the function with an adapter with broyden result, iterations = broyden2(fixed_point_1_guess, fixed_point_to_residual(fixed_point_1_func), jac=None, skip_J=True, xtol=1e-9, maxiter=1000) assert iterations < 50 # 38 last time assert_close1d(result, fixed_point_1_expect, rtol=1e-9) # Test it with standard newton forms: def basic_system(inputs): x, y = inputs return [x**2 + y**2 - 4, exp(x) - y] # Initial guess, note it does not converge with damping of 1. It seems many functions will not converge no matter the function. x0 = [1.0, 1.0] solution, iterations = fixed_point( f=residual_to_fixed_point(basic_system), x0=x0, xtol=1e-10, maxiter=1000, damping=.3, ) assert iterations < 80 # 74 last check assert_close1d(basic_system(solution), [0, 0], atol=1e-6) # Check the other methods converge solution, iterations = fixed_point_aitken( f=residual_to_fixed_point(basic_system), x0=x0, xtol=1e-10, maxiter=1000, damping=0.3, acc_damping=0.4,# 44 last run iterations ) assert iterations < 100 assert_close1d(basic_system(solution), [0, 0], atol=1e-6) # Check the other methods converge solution, iterations = fixed_point_gdem( f=residual_to_fixed_point(basic_system), x0=x0, xtol=1e-10, maxiter=1000, damping=0.3, acc_damping=0.3, # 81 last run ) assert iterations < 100 assert_close1d(basic_system(solution), [0, 0], atol=1e-6) # # Anderson acceleration does something different, always solves for the same set of residuals solution, iterations = fixed_point_anderson( f=(basic_system), x0=x0, xtol=1e-6, maxiter=1000, window_size=2, reg=0, mixing_param=0.5, damping=0.3, initial_iterations=0 ) assert iterations < 15 # 25 last time # assert_close1d(basic_system(solution), [0, 0], atol=1e-6) # def to_fixed_point(inputs): # x, y = inputs # new_y = exp(x) # new_x = (abs(4 - y*y))**0.5 # return [new_x, new_y] # # Initial guess # x0 = [1.0, 1] # solution, iterations = fixed_point_anderson( # f=(to_fixed_point), # x0=x0, # xtol=1e-6, # maxiter=1000, # window_size=2, reg=0, mixing_param=0.5, # damping=0.3, # initial_iterations=0 # ) # assert iterations < 15 # 25 last time # # assert_close1d(basic_system(solution), [0, 0], atol=1e-6) def test_SolverInterface_fixed_point(): # Largely from Yoel's Flexsolve which is excellent, and wikipedia sample code https://en.wikipedia.org/wiki/Aitken%27s_delta-squared_process # Works really nice, few parameters not exposed still though solver = SolverInterface(method='fixed_point', objf=fixed_point_to_residual(fixed_point_1_func), maxiter=1000, xtol=1e-9) ans = solver.solve(fixed_point_1_guess) assert_close1d(ans, fixed_point_1_expect, rtol=1e-6) solver = SolverInterface(method='fixed_point_aitken', objf=fixed_point_to_residual(fixed_point_1_func), maxiter=1000, xtol=1e-12) ans = solver.solve(fixed_point_1_guess) assert solver.fval_iter == 79 assert_close1d(ans, fixed_point_1_expect, rtol=1e-6) solver = SolverInterface(method='fixed_point_gdem', objf=fixed_point_to_residual(fixed_point_1_func), maxiter=1000, xtol=1e-12) ans = solver.solve(fixed_point_1_guess) assert solver.fval_iter == 37 assert_close1d(ans, fixed_point_1_expect, rtol=1e-12) def test_isclose(): assert not isclose(1.0, 2.0, rel_tol=0.0, abs_tol=0.0) assert not isclose(1.0, 1+1e-14, rel_tol=0.0, abs_tol=0.0) assert isclose(1.0, 1+1e-14, rel_tol=1e-13, abs_tol=0.0) assert isclose(1.0, 1+1e-14, rel_tol=0.0, abs_tol=1e-14) assert not isclose(1.0, 1+1e-14, rel_tol=0.0, abs_tol=9e-15) assert isclose(1e300, 1e300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) assert isclose(1e-300, 1e-300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) assert isclose(-1e300, -1e300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) assert isclose(-1e-300, -1e-300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) assert not isclose(1.0, 2.0, rel_tol=0.0, abs_tol=0.0) assert not isclose(1.0, 1+1e-14, rel_tol=0.0, abs_tol=0.0) assert isclose(1.0, 1+1e-14, rel_tol=1e-13, abs_tol=0.0) assert isclose(1.0, 1+1e-14, rel_tol=0.0, abs_tol=1e-14) assert not isclose(1.0, 1+1e-14, rel_tol=0.0, abs_tol=9e-15) assert isclose(1e300, 1e300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) assert isclose(1e-300, 1e-300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) assert isclose(-1e300, -1e300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) assert isclose(-1e-300, -1e-300*(1+1e-14), rel_tol=2e-14, abs_tol=0.0) def test_trunc_exp_numpy(): dat = np.array([-1000, 1e4, -1.2, 0.0, 1.0, 1000.0]) expect = [0.0, 1.7976931348622732e+308, 0.30119421191220214, 1.0, 2.718281828459045, 1.7976931348622732e+308] calc = trunc_exp_numpy(dat) assert_close1d(calc, expect, rtol=1e-13) assert not np.any(np.isnan(expect)) assert not np.any(np.isinf(expect)) out_array = np.zeros(6) calc = trunc_exp_numpy(dat, out=out_array) assert_close1d(calc, expect, rtol=1e-13) assert not np.any(np.isnan(expect)) assert not np.any(np.isinf(expect)) assert out_array is calc def test_trunc_log_numpy(): dat = np.array([0.0, 1e4, 0.0, 1.0, 1000.0]) expect = [-744.4400719213812, 9.210340371976184, -744.4400719213812, 0.0, 6.907755278982137] calc = trunc_log_numpy(dat) assert_close1d(calc, expect, rtol=1e-13) assert not np.any(np.isnan(expect)) assert not np.any(np.isinf(expect)) out_array = np.zeros(5) calc = trunc_log_numpy(dat, out=out_array) assert_close1d(calc, expect, rtol=1e-13) assert not np.any(np.isnan(expect)) assert not np.any(np.isinf(expect)) assert out_array is calc dat = np.array([[0.0, 1e4, 1.0], [2.0, 0.0, 1000.0]]) expect = [[-744.4400719213812, 9.210340371976184, 0.0], [0.6931471805599453, -744.4400719213812, 6.907755278982137]] calc = trunc_log_numpy(dat) assert_close2d(calc, expect, rtol=1e-13) assert not np.any(np.isnan(expect)) assert not np.any(np.isinf(expect)) out_array = np.zeros(6).reshape((2, 3)) calc = trunc_log_numpy(dat, out=out_array) assert_close2d(calc, expect, rtol=1e-13) assert not np.any(np.isnan(expect)) assert not np.any(np.isinf(expect)) assert out_array is calc def test_transpose(): from fluids.numerics import transpose l=[[1,2,3],[4,5,6],[7,8,9]] out = transpose(l) assert_close2d(out, [[1, 4, 7], [2, 5, 8], [3, 6, 9]]) def test_is_increasing(): assert is_increasing([1, 2, 3, 4]) assert not is_increasing([1, 1, 2, 3]) assert not is_increasing([4, 3, 2, 1]) assert is_increasing([1000000, 1000001, 1000002, 1000003]) assert not is_increasing([1000003, 1000002, 1000001, 1000000]) assert is_increasing([1.1, 2.2, 3.3, 4.4]) assert not is_increasing([4.4, 3.3, 2.2, 1.1]) assert is_increasing([1, 2.5, 3, 4.5]) assert not is_increasing([4.5, 3, 2.5, 1]) assert is_increasing([1000000]) assert is_increasing([3.14159]) assert is_increasing([-4, -3, -2, -1]) assert not is_increasing([-1, -2, -3, -4]) assert is_increasing([]) assert is_increasing([1, 2]) assert not is_increasing([2, 1]) def test_py_lambertw(): from fluids.numerics import py_lambertw # from scipy.special import lambertw as py_lambertw """ Test the basic functionality of the Lambert W function """ # Check a simple positive value assert abs(py_lambertw(0).real - 0) < 0.22 # Check a value greater than 1 result = py_lambertw(exp(1) - 1).real # LambertW(e-1) should be 1 assert_close(abs(result - 1), 0.21570504021558445) # Test the boundary at -1/e minus_one_over_e = -1/exp(1) res = abs(py_lambertw(minus_one_over_e, k=-1).real + 1) def test_hessian(): def f(x): return x[0]**2 + x[1]**2 x0 = [1.0, 2.0] result = hessian(f, x0, perturbation=1e-4) expected = [[2.0, 0.0], [0.0, 2.0]] assert_close2d(result, expected) def f(x): return 1000*sin(x[0]) + 4000*cos(x[1]) x0 = [.3, .7] result = hessian(f, x0, perturbation=1e-5) expected = [[-296.1312838005886, 0.0], [0.0, -3055.5183310366906]] assert_close2d(result, expected, rtol=0.1) def f(x): return x[0]**2 + 2*x[1]**2 + x[0]*x[1] x0 = [1.0, 1.0] result = hessian(f, x0, full=False, perturbation=1e-4) expected = [[2.0], [1.0, 4.0]] assert_close1d(result[0], expected[0]) assert_close1d(result[1], expected[1]) def f(x): return [x[0]**2 - x[1], x[0]*x[1]] x0 = [1.0, 1.0] result = hessian(f, x0, scalar=False, perturbation=1e-4) expected = [[[2.0, 0.0], [0.0, 1.0]], [[0.0, 1.0], [0.0, 0.0]]] assert_close3d(result, expected, rtol=1e-5)