'''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 log10 import pytest import fluids from fluids.numerics import assert_close, assert_close1d, assert_close2d, isclose, logspace from fluids.particle_size_distribution import ( ISO_3310_1_R10, ISO_3310_1_R20, ISO_3310_1_R20_3, ISO_3310_1_R40_3, ASTM_E11_sieves, ISO_3310_1_sieves, ParticleSizeDistribution, PSDCustom, PSDInterpolated, PSDLognormal, cdf_Gates_Gaudin_Schuhman, cdf_lognormal, cdf_Rosin_Rammler, pdf_Gates_Gaudin_Schuhman, pdf_Gates_Gaudin_Schuhman_basis_integral, pdf_lognormal, pdf_lognormal_basis_integral, pdf_Rosin_Rammler, pdf_Rosin_Rammler_basis_integral, psd_spacing, ) try: from random import uniform except: pass def test_ASTM_E11_sieves(): sieves = ASTM_E11_sieves.values() tot = sum([i.d_wire for i in sieves]) assert_close(tot, 0.105963384) tot = sum([i.opening for i in sieves]) assert_close(tot, 0.9876439999999999) assert len(ASTM_E11_sieves) == 56 # Test but do not validate these properties tot = 0.0 for attr in ['Y_variation_avg', 'X_variation_max', 'max_opening', 'd_wire', 'd_wire_min', 'd_wire_max', 'opening', 'opening_inch']: tot += sum(getattr(i, attr) for i in sieves) def test_ISO_3310_2_sieves(): sieves = ISO_3310_1_sieves.values() tot = sum([i.d_wire for i in sieves]) assert_close(tot, 0.17564599999999997) tot = sum([i.opening for i in sieves]) assert_close(tot, 1.5205579999999994) assert len(ISO_3310_1_sieves) == 99 # Test but do not validate these properties tot = 0.0 for attr in ['Y_variation_avg', 'X_variation_max', 'd_wire', 'd_wire_min', 'd_wire_max', 'opening']: tot += sum(getattr(i, attr) for i in sieves) for l in [ISO_3310_1_R20_3, ISO_3310_1_R20, ISO_3310_1_R10, ISO_3310_1_R40_3]: for i in l: assert i.designation in ISO_3310_1_sieves def test_ParticleSizeDistribution_basic(): ds = [240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532] numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1] number_fractions = [0.010640039286298903, 0.01947945653953184, 0.03797675560648224, 0.06711409395973154, 0.102962841708954, 0.13897528237027337, 0.16205598297593715, 0.160582746767065, 0.13504665247994763, 0.09477819610410869, 0.048616794892781146, 0.01816991324275659, 0.0034375511540350305, 0.0001636929120969062] length_fractions = [0.0022265080273913248, 0.005405749400984079, 0.013173675010801534, 0.02909808308708846, 0.05576732372469186, 0.09403390879219536, 0.1370246122004729, 0.16966553692650058, 0.17831420382670332, 0.15641421494054603, 0.10028800800464328, 0.046849963047687335, 0.011078803825079166, 0.0006594091852147985] area_fractions = [0.0003643458522227456, 0.0011833425086503686, 0.0036047198267710797, 0.009951607879295004, 0.023826910138492176, 0.05018962198499494, 0.09139246506396961, 0.1414069073893575, 0.18572285033413602, 0.20362023102799823, 0.16318760564859225, 0.09528884410476045, 0.028165197280747324, 0.0020953509600122053] fractions = [4.8560356399310335e-05, 0.00021291794698947167, 0.0008107432330218852, 0.0027975134942445257, 0.00836789808490677, 0.02201901107895143, 0.05010399231412809, 0.0968727835386488, 0.15899879607747244, 0.2178784903712532, 0.21825921197532888, 0.159302671180342, 0.05885464261922434, 0.0054727677290887945] fraction_cdf = [4.856035639931034e-05, 0.0002614783033887821, 0.0010722215364106676, 0.003869735030655194, 0.012237633115561966, 0.0342566441945134, 0.0843606365086415, 0.18123342004729032, 0.34023221612476284, 0.5581107064960161, 0.7763699184713451, 0.9356725896516871, 0.9945272322709114, 1.0000000000000002] area_cdf = [0.00036434585222274563, 0.0015476883608731143, 0.005152408187644195, 0.015104016066939202, 0.038930926205431385, 0.08912054819042634, 0.18051301325439598, 0.3219199206437535, 0.5076427709778896, 0.7112630020058879, 0.8744506076544801, 0.9697394517592406, 0.9979046490399879, 1.0] length_cdf = [0.0022265080273913248, 0.007632257428375404, 0.020805932439176937, 0.0499040155262654, 0.10567133925095726, 0.1997052480431526, 0.3367298602436255, 0.506395397170126, 0.6847096009968294, 0.8411238159373755, 0.9414118239420188, 0.9882617869897061, 0.9993405908147853, 1.0000000000000002] number_cdf = [0.010640039286298902, 0.030119495825830737, 0.06809625143231296, 0.13521034539204452, 0.2381731871009985, 0.3771484694712718, 0.5392044524472088, 0.6997871992142738, 0.8348338516942214, 0.9296120477983301, 0.9782288426911112, 0.9963987559338677, 0.9998363070879027, 0.9999999999999997] opts = [ {'fractions': numbers, 'cdf': False, 'order': 0}, {'fractions': number_fractions, 'cdf': False, 'order': 0}, {'fractions': length_fractions, 'cdf': False, 'order': 1}, {'fractions': area_fractions, 'cdf': False, 'order': 2}, {'fractions': fractions, 'cdf': False, 'order': 3}, {'fractions': fraction_cdf, 'cdf': True, 'order': 3}, {'fractions': area_cdf, 'cdf': True, 'order': 2}, {'fractions': length_cdf, 'cdf': True, 'order': 1}, {'fractions': number_cdf, 'cdf': True, 'order': 0}] for opt in opts: asme_e799 = ParticleSizeDistribution(ds=ds, **opt) d10 = asme_e799.mean_size(1, 0) assert_close(d10, 1459.3725650679328) d21 = asme_e799.mean_size(2, 1) assert_close(d21, 1857.7888572055529) d20 = asme_e799.mean_size(2, 0) assert_close(d20, 1646.5740462835831) d32 = asme_e799.mean_size(3, 2) assert_close(d32, 2269.3210317450453) # This one is rounded to 2280 in ASME - weird d31 = asme_e799.mean_size(3, 1) assert_close(d31, 2053.2703977309357) # This one is rounded to 2060 in ASME - weird d30 = asme_e799.mean_size(3, 0) assert_close(d30, 1832.39665294744) d43 = asme_e799.mean_size(4, 3) assert_close(d43, 2670.751954612969) # The others are, rounded to the nearest 10, correct. # There's something weird about the end points of some intermediate values of # D3 and D4. Likely just rounding issues. vol_percents_exp = [0.005, 0.021, 0.081, 0.280, 0.837, 2.202, 5.010, 9.687, 15.900, 21.788, 21.826, 15.930, 5.885, 0.547] assert vol_percents_exp == [round(i*100, 3) for i in asme_e799.fractions] assert_close1d(asme_e799.fractions, fractions) assert_close1d(asme_e799.number_fractions, number_fractions) # i, i distributions d00 = asme_e799.mean_size(0, 0) assert_close(d00, 1278.7057976023061) d11 = asme_e799.mean_size(1, 1) assert_close(d11, 1654.6665309027303) d22 = asme_e799.mean_size(2, 2) assert_close(d22, 2054.3809583432208) d33 = asme_e799.mean_size(3, 3) assert_close(d33, 2450.886241250387) d44 = asme_e799.mean_size(4, 4) assert_close(d44, 2826.0471682278476) vssa = asme_e799.vssa assert_close(vssa, 0.0026656187302839165) def test_pdf_lognormal(): import scipy.stats pdf = pdf_lognormal(d=1E-4, d_characteristic=1E-5, s=1.1) assert_close(pdf, 405.5420921156425, rtol=1E-12) pdf_sp = scipy.stats.lognorm.pdf(x=1E-4/1E-5, s=1.1)/1E-5 assert_close(pdf_sp, pdf) assert 0.0 == pdf_lognormal(d=0, d_characteristic=1E-5, s=1.1) # Check we can get down almost to zero pdf = pdf_lognormal(d=3.7E-24, d_characteristic=1E-5, s=1.1) assert_close(pdf, 4.842842147909424e-301, atol=1e-200) def test_cdf_lognormal(): import scipy.stats cdf = cdf_lognormal(d=1E-4, d_characteristic=1E-5, s=1.1) assert_close(cdf, 0.98183698757981763) cdf_sp = scipy.stats.lognorm.cdf(x=1E-4/1E-5, s=1.1) assert_close(cdf, cdf_sp) assert cdf_lognormal(d=1e300, d_characteristic=1E-5, s=1.1) == 1.0 assert cdf_lognormal(d=0, d_characteristic=1E-5, s=1.1) == 0.0 def test_pdf_lognormal_basis_integral(): ans = pdf_lognormal_basis_integral(d=1E-4, d_characteristic=1E-5, s=1.1, n=-2) assert_close(ans, 56228306549.263626) # Some things: ans = pdf_lognormal_basis_integral(d=1E-100, d_characteristic=1E-5, s=1.1, n=-2) ans2 = pdf_lognormal_basis_integral(d=1E-120, d_characteristic=1E-5, s=1.1, n=-2) assert_close(ans, ans2, rtol=1E-12) # Couldn't get the limit for pdf_lognormal_basis_integral when d = 0 # # with Sympy Or wolfram @pytest.mark.scipy @pytest.mark.fuzz @pytest.mark.slow def test_pdf_lognormal_basis_integral_fuzz(): from scipy.integrate import quad # The actual integral testing analytical_vales = [] numerical_values = [] for n in [-3, -2, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]: # Errors at -2 (well, prevision loss anyway) for d_max in [1e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 6e-2, 7e-2, 8e-2, 1e-1]: d_max = d_max/100 # Make d smaller analytical = (pdf_lognormal_basis_integral(d_max, d_characteristic=1E-5, s=1.1, n=n) - pdf_lognormal_basis_integral(1e-20, d_characteristic=1E-5, s=1.1, n=n)) to_int = lambda d : d**n*pdf_lognormal(d, d_characteristic=1E-5, s=1.1) points = logspace(log10(d_max/1000), log10(d_max*.999), 40) numerical = quad(to_int, 1e-9, d_max, points=points)[0] # points=points analytical_vales.append(analytical) numerical_values.append(numerical) assert_close1d(analytical_vales, numerical_values, rtol=2E-6) def test_cdf_Gates_Gaudin_Schuhman(): ''' d, d_max, n, m = symbols('d, d_max, n, m') expr = (d/d_max)**n pdf = diff(expr, d) ''' cdf = cdf_Gates_Gaudin_Schuhman(d=2E-4, d_characteristic=1E-3, m=2.3) assert_close(cdf, 0.024681354508800397) cdf = cdf_Gates_Gaudin_Schuhman(d=1.01e-3, d_characteristic=1E-3, m=2.3) assert_close(cdf, 1) def test_pdf_Gates_Gaudin_Schuhman(): pdf = pdf_Gates_Gaudin_Schuhman(d=2E-4, d_characteristic=1E-3, m=2.3) assert_close(pdf, 283.8355768512045) pdf = pdf_Gates_Gaudin_Schuhman(d=2E-3, d_characteristic=1E-3, m=2.3) assert_close(pdf, 0) def test_pdf_Gates_Gaudin_Schuhman_basis_integral(): ans = pdf_Gates_Gaudin_Schuhman_basis_integral(d=0, d_characteristic=1e-3, m=2.3, n=5) assert_close(ans, 0) ans = pdf_Gates_Gaudin_Schuhman_basis_integral(d=1e-3, d_characteristic=1e-3, m=2.3, n=5) assert_close(ans, 3.1506849315068495e-16) @pytest.mark.scipy @pytest.mark.slow @pytest.mark.fuzz def test_pdf_Gates_Gaudin_Schuhman_basis_integral_fuzz(): from scipy.integrate import quad """Note: Test takes 10x longer with the quad/points. """ analytical_vales = [] numerical_values = [] for n in [-1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]: # Errors at -2 (well, prevision loss anyway) for d_max in [1e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 6e-2, 7e-2, 8e-2, 1e-1]: d_max = d_max/100 # d cannot be larger than d_max analytical = (pdf_Gates_Gaudin_Schuhman_basis_integral(d_max, 1E-3, 2.3, n) - pdf_Gates_Gaudin_Schuhman_basis_integral(1E-20, 1E-3, 2.3, n)) to_int = lambda d : d**n*pdf_Gates_Gaudin_Schuhman(d, 1E-3, 2.3) numerical = quad(to_int, 1E-20, d_max)[0] # points=points analytical_vales.append(analytical) numerical_values.append(numerical) # The precision here is amazing actually, 1e-14 passes # assert_close1d(analytical, numerical, rtol=1E-7) assert_close1d(analytical_vales, numerical_values, rtol=1E-7) def test_cdf_Rosin_Rammler(): cdf = cdf_Rosin_Rammler(5E-2, 200.0, 2.0) assert_close(cdf, 0.3934693402873667) def test_pdf_Rosin_Rammler(): ''' from sympy import * d, k, n = symbols('d, k, n') model = 1 - exp(-k*d**n) print(latex(diff(model, d))) ''' from scipy.integrate import quad pdf = pdf_Rosin_Rammler(1E-3, 200.0, 2.0) assert_close(pdf, 0.3999200079994667) # quad to_quad = lambda d: pdf_Rosin_Rammler(d, 200.0, 2.0) cdf_int = quad(to_quad, 0, 5e-2)[0] cdf_known = cdf_Rosin_Rammler(5E-2, 200, 2) assert_close(cdf_int, cdf_known) assert_close(1, quad(to_quad, 0, 5)[0]) assert 0 == pdf_Rosin_Rammler(0, 200, 2) def test_pdf_Rosin_Rammler_basis_integral(): ans = pdf_Rosin_Rammler_basis_integral(5E-2, 200.0, 2.0, 3) assert_close(ans, -0.00045239898439007338) # Test no error pdf_Rosin_Rammler_basis_integral(0, 200, 2, 3) @pytest.mark.slow @pytest.mark.fuzz def test_pdf_Rosin_Rammler_basis_integral_fuzz(): from scipy.integrate import quad for n in [1.0, 2.0, 3.0]: # Lower d_maxes have for d_max in [ 1e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 6e-2, 7e-2, 8e-2, 1e-1]: analytical = (pdf_Rosin_Rammler_basis_integral(d_max, 200, 2, n) - pdf_Rosin_Rammler_basis_integral(1E-20, 200, 2, n)) to_int = lambda d : d**n*pdf_Rosin_Rammler(d, 200, 2) points = logspace(log10(d_max/2000), log10(d_max*.999), 30) numerical = quad(to_int, 1E-20, d_max, points=points)[0] assert_close(analytical, numerical, rtol=1E-5) def testPSDLognormal_meshes(): a = PSDLognormal(s=0.5, d_characteristic=5E-6) ds_expect = [5.011872336272722e-07, 6.309573444801932e-07, 7.943282347242815e-07, 1e-06] ds = a.ds_discrete(d_max=1e-6, method='R10', pts=4) assert_close1d(ds_expect, ds) ds = a.ds_discrete(d_min=1e-6, method='R10', pts=4) assert_close1d(ds, [1e-06, 1.2589254117941672e-06, 1.5848931924611134e-06, 1.9952623149688796e-06]) @pytest.mark.fuzz @pytest.mark.slow def test_PSDLognormal_mean_sizes_numerical(): '''Takes like 1 second. Should not run normally. Not how things should be done, just a proof of concept. ''' # ISO standard example, done numerically a = PSDLognormal(s=0.5, d_characteristic=5E-6) ds = a.ds_discrete(d_max=1, pts=1E5) fractions = a.fractions_discrete(ds) disc = ParticleSizeDistribution(ds=ds, fractions=fractions, order=3) d20 = disc.mean_size(2, 0) assert_close(d20, 3.033E-6, rtol=0, atol=1E-9) d10 = disc.mean_size(1, 0) assert_close(d10, 2.676E-6, rtol=0, atol=1E-9) d21 = disc.mean_size(2, 1) assert_close(d21, 3.436E-6, rtol=0, atol=1E-9) d32 = disc.mean_size(3, 2) # Does match, need 6E6 pts to get the last digit to round right assert_close(d32, 4.412E-6, rtol=0, atol=1E-9) d43 = disc.mean_size(4, 3) assert_close(d43, 5.666E-6, rtol=0, atol=1E-9) d33 = disc.mean_size(3.0, 3.0) assert_close(d33, 5.000E-6, rtol=0, atol=1E-9) d00 = disc.mean_size(0.0, 0.0) assert_close(d00, 2.362E-6, rtol=0, atol=1E-9) @pytest.mark.skipif(fluids.numerics.IS_PYPY, reason="numerical issues with adaptive integrator") # Custom integrator can't handle this @pytest.mark.slow def test_PSDCustom_mean_sizes_numerical(): import scipy.stats distribution = scipy.stats.lognorm(s=0.5, scale=5E-6) disc = PSDCustom(distribution) d20 = disc.mean_size(2, 0) assert_close(d20, 3.0326532985631672e-06, rtol=5e-7) d10 = disc.mean_size(1, 0) assert_close(d10, 2.6763071425949508e-06, rtol=1e-6) d21 = disc.mean_size(2, 1) assert_close(d21, 3.4364463939548618e-06, rtol=1e-7) d32 = disc.mean_size(3, 2) assert_close(d32, 4.4124845129229773e-06, rtol=1E-7) d43 = disc.mean_size(4, 3) assert_close(d43, 5.6657422653341318e-06, rtol=1E-3) def test_PSDLognormal_mean_sizes_analytical(): disc = PSDLognormal(s=0.5, d_characteristic=5E-6) d20 = disc.mean_size(2, 0) assert_close(d20, 3.033E-6, rtol=0, atol=1E-9) assert_close(d20, 3.0326532985631672e-06, rtol=1E-12) assert_close(d20, disc.mean_size_ISO(2, 0), rtol=1E-12) d10 = disc.mean_size(1, 0) assert_close(d10, 2.676E-6, rtol=0, atol=1E-9) assert_close(d10, 2.6763071425949508e-06, rtol=1E-12) assert_close(d10, disc.mean_size_ISO(1, 0), rtol=1E-12) d21 = disc.mean_size(2, 1) assert_close(d21, 3.436E-6, rtol=0, atol=1E-9) assert_close(d21, 3.4364463939548618e-06, rtol=1E-12) assert_close(d21, disc.mean_size_ISO(1, 1), rtol=1E-12) d32 = disc.mean_size(3, 2) assert_close(d32, 4.412E-6, rtol=0, atol=1E-9) assert_close(d32, 4.4124845129229773e-06, rtol=1E-12) assert_close(d32, disc.mean_size_ISO(1, 2), rtol=1E-12) d43 = disc.mean_size(4, 3) assert_close(d43, 5.666E-6, rtol=0, atol=1E-9) assert_close(d43, 5.6657422653341318e-06, rtol=1E-12) assert_close(d43, disc.mean_size_ISO(1, 3), rtol=1E-12) assert_close(disc.vssa, 1359778.1436801916) # There guys - need more work d33 = disc.mean_size(3.0, 3.0) assert_close(d33, 5.000E-6, rtol=0, atol=1E-6) # d00 = disc.mean_size(0.0, 0.0) assert_close(d00, 2.362E-6, rtol=1e-4) def test_PSDLognormal_ds_discrete(): ds_discrete_expect = [2.4920844782646124e-07, 5.870554386661556e-07, 1.3829149496067687e-06, 3.2577055451375563e-06, 7.674112874286059e-06, 1.8077756749742233e-05, 4.258541598963051e-05, 0.0001003176268004454] dist = PSDLognormal(s=0.5, d_characteristic=5E-6) # Test the defaults assert_close1d(dist.ds_discrete(pts=8), ds_discrete_expect, rtol=1e-4) ds_discrete_expect = [1e-07, 1.389495494373139e-07, 1.9306977288832497e-07, 2.6826957952797275e-07, 3.727593720314938e-07, 5.179474679231213e-07, 7.196856730011514e-07, 1e-06] # Test end and minimum points assert_close1d(dist.ds_discrete(d_min=1e-7, d_max=1e-6, pts=8), ds_discrete_expect, rtol=1e-12) def test_PSDLognormal_ds_discrete(): # Test the cdf discrete dist = PSDLognormal(s=0.5, d_characteristic=5E-6) ds = dist.ds_discrete(d_min=1e-7, d_max=1e-6, pts=8) ans = dist.fractions_discrete(ds) fractions_expect = [2.55351295663786e-15, 3.831379657981415e-13, 3.762157252396037e-11, 2.41392961175535e-09, 1.01281244724305e-07, 2.7813750147487326e-06, 5.004382447515443e-05, 0.00059054208024234] assert_close1d(fractions_expect, ans, rtol=1e-3) def test_PSDLognormal_dn(): disc = PSDLognormal(s=0.5, d_characteristic=5E-6) # Test input of 1 ans = disc.dn(1) # The answer can vary quite a lot near the end, so it is safest just to # compare with the reverse, plugging it back to cdf assert_close(disc.cdf(ans), 1, rtol=1E-12) # assert_close(ans, 0.0002964902595794474) # Test zero input assert_close(disc.dn(0), 0) # Test 50% input ans = disc.dn(.5) assert_close(ans, 5.0e-06, rtol=1E-6) with pytest.raises(Exception): disc.dn(1.5) with pytest.raises(Exception): disc.dn(-.5) # Other orders of n - there is no comparison data for this yet!! assert_close(disc.pdf(1E-5), disc.pdf(1E-5, 3)) assert_close(disc.pdf(1E-5, 2), 13468.122877854335) assert_close(disc.pdf(1E-5, 1), 4628.2482296943508) assert_close(disc.pdf(1E-5, 0), 1238.6613794833427) # Some really large s tests - found some issues with this dist = PSDLognormal(s=4, d_characteristic=5E-6) assert_close(dist.dn(1e-15), 8.220922763476676e-20, rtol=1e-1) assert_close(dist.dn(.99999999), 28055.285560763594) assert_close(dist.dn(1e-9), 1.904197766691136e-16, rtol=1e-4) assert_close(dist.dn(1-1e-9), 131288.88851649483, rtol=1e-4) def test_PSDLognormal_dn_order_0_1_2(): '''Simple point to test where the order of n should be 0 Yes, the integrals need this many points (which makes them slow) to get the right accuracy. They've been tested and reduced already quite a bit. ''' from scipy.integrate import quad # test 2, 0 -> 2, 0 disc = PSDLognormal(s=0.5, d_characteristic=5E-6) to_int = lambda d: d**2*disc.pdf(d=d, n=0) points = [5E-6*i for i in logspace(log10(.1), log10(50), 8)] ans_numerical = (quad(to_int, 1E-7, 5E-3, points=points)[0])**0.5 ans_analytical = 3.0326532985631672e-06 # The integral is able to give over to decimals! assert_close(ans_numerical, ans_analytical, rtol=1E-10) # test 2, 1 -> 1, 1 integrated pdf to_int = lambda d: d*disc.pdf(d=d, n=1) ans_numerical = (quad(to_int, 1E-7, 5E-3, points=points)[0])**1 ans_analytical = 3.4364463939548618e-06 assert_close(ans_numerical, ans_analytical, rtol=1E-10) # test 3, 2 -> 1, 2 integrated pdf to_int = lambda d: d*disc.pdf(d=d, n=2) ans_numerical = (quad(to_int, 1E-7, 5E-3, points=points)[0])**1 ans_analytical = 4.4124845129229773e-06 assert_close(ans_numerical, ans_analytical, rtol=1E-8) def test_PSDLognormal_cdf_orders(): # Test cdf of different orders a bunch disc = PSDLognormal(s=0.5, d_characteristic=5E-6) # 16 x 4 = 64 points # had 1e-7 here too as a diameter but too many numerical issues, too sensitive to rounding # errors ds = [2E-6, 3E-6, 4E-6, 5E-6, 6E-6, 7E-6, 1E-5, 2E-5, 3E-5, 5E-5, 7E-5, 1E-4, 2E-4, 4E-4, 1E-3] ans_expect = [[ 0.36972511868508068, 0.68379899882263917, 0.8539928088656249, 0.93319279873114203, 0.96888427729983861, 0.98510775165387254, 0.99805096305713792, 0.9999903391682271, 0.99999981474719135, 0.99999999948654394, 0.99999999999391231, 0.99999999999996592, 1.0, 1.0, 1.0], [ 0.20254040832522924, 0.49136307673913149, 0.71011232639847854, 0.84134474606854293, 0.91381737643345484, 0.95283088619207579, 0.99149043874391107, 0.99991921875653167, 0.99999771392273817, 0.99999998959747816, 0.99999999982864851, 0.99999999999863987, 1.0, 1.0, 1.0000000000000002], [ 0.091334595732478097, 0.30095658738958564, 0.52141804648990697, 0.69146246127401301, 0.80638264936531323, 0.87959096325267294, 0.9703723506333426, 0.999467162897961, 0.99997782059383122, 0.99999983475152954, 0.99999999622288382, 0.99999999995749711, 0.99999999999999833, 1.0000000000000002, 1.0000000000000002], [ 0.033432418408916864, 0.15347299656473007, 0.3276949357115424, 0.5, 0.64231108623683952, 0.74950869138681098, 0.91717148099830148, 0.99721938213769046, 0.99983050191355338, 0.99999793935660408, 0.99999993474010451, 0.99999999896020164, 0.99999999999991951, 1.0, 1.0]] calc = [] for n in range(0, 4): calc.append([disc.cdf(i, n=n) for i in ds]) assert_close2d(ans_expect, calc, rtol=1E-9) def test_PSDLognormal_cdf_vs_pdf(): from scipy.integrate import quad # test PDF against CDF disc = PSDLognormal(s=0.5, d_characteristic=5E-6) ans_calc = [] ans_expect = [] for i in range(5): # Pick a random start start = uniform(0, 1) end = uniform(start, 1) d_start = disc.dn(start) d_end = disc.dn(end) delta = disc.cdf(d_end) - disc.cdf(d_start) delta_numerical = quad(disc.pdf, d_start, d_end)[0] ans_calc.append(delta_numerical) ans_expect.append(delta) assert_close1d(ans_calc, ans_expect) def test_PSD_lognormal_truncated(): psd = PSDLognormal(s=0.5, d_characteristic=5E-6, d_max=1e-5, d_min=1e-6, order=3) assert_close(psd.dn(1), 1e-5) assert_close(psd.dn(0), 1e-6) assert_close(psd.cdf(psd.dn(0.5)), 0.5) # No longer at d_characteristic # Check the cdf limits of the truncated distribution ds = [1e-8, 1e-7, 1e-6, 2e-6, 5e-6, 9e-6, 1e-5, 1e-4, 1e-3] ans = psd.fractions_discrete(ds) ans_expect = [0.0, 0.0, 0.0, 0.03577517221380477, 0.5090598176028703, 0.41473614196545805, 0.04042886821786684, 0.0, 0.0] assert sum(ans) == 1 assert_close1d(ans, ans_expect) pdfs = [psd.pdf(i) for i in ds] pdfs_expect = [0.0, 0.0, 4896.6230234795203, 81190.803665720552, 174110.24040947447, 48468.576095204902, 33302.599459012476, 0.0, 0.0] assert_close1d(pdfs, pdfs_expect) # Check the mean_size calculations - this is right assert_close(psd.mean_size(3, 2), 4.1817574249337556e-06) @pytest.mark.fuzz @pytest.mark.slow def test_PSD_lognormal_truncated_mean_size(): from scipy.integrate import quad psd = PSDLognormal(s=0.5, d_characteristic=5E-6, d_max=1e-5, d_min=1e-6, order=3) def psd_mean_size_numeric_truncated(psd, p, q): pow1 = q - psd.order to_int = lambda d : d**pow1*psd.pdf(d) denominator = quad(to_int, psd.d_min, psd.d_max)[0] # denominator = self._pdf_basis_integral_definite(d_min=self.d_minimum, d_max=self.d_excessive, n=pow1) root_power = p - q pow3 = p - psd.order # numerator = self._pdf_basis_integral_definite(d_min=self.d_minimum, d_max=self.d_excessive, n=pow3) to_int = lambda d : d**pow3*psd.pdf(d) numerator = quad(to_int, psd.d_min, psd.d_max)[0] return (numerator/denominator)**(1.0/(root_power)) # It really works! for p in range(-3, 4): for q in range(-3, 4): if p != q: assert_close(psd.mean_size(p, q), psd_mean_size_numeric_truncated(psd, p, q)) @pytest.mark.slow def test_PSD_PSDlognormal_area_length_count(): '''Compare the average difference between the analytical values for a lognormal distribution with those of a discretized form of it.Note simply adding more points did not tend to help reduce the error. For the particle count case, 700 points has the lowest error. fractions_discrete is still the slowest part. ''' import numpy as np dist = PSDLognormal(s=0.5, d_characteristic=5E-6) ds = dist.ds_discrete(pts=700, d_min=2E-7, d_max=1E-4) fractions = dist.fractions_discrete(ds) psd = ParticleSizeDistribution(ds=ds, fractions=fractions, order=3) # Trim a few at the start and end ans = np.array(psd.number_fractions)[5:-5]/np.array(dist.fractions_discrete(ds, n=0))[5:-5] avg_err = sum(np.abs(ans - 1.0))/len(ans) assert 5E-4 > avg_err ans = np.array(psd.length_fractions)[5:-5]/np.array(dist.fractions_discrete(ds, n=1))[5:-5] avg_err = sum(np.abs(ans - 1.0))/len(ans) assert 1E-4 > avg_err ans = np.array(psd.area_fractions)[5:-5]/np.array(dist.fractions_discrete(ds, n=2))[5:-5] avg_err = sum(np.abs(ans - 1.0))/len(ans) assert 1E-4 > avg_err def test_PSDInterpolated_pchip(): '''For this test, ds is the same length as fractions, and we begin the series with the zero point. Half the test is spend on the `dn` solver tests, and the other half is just that these tests are slow. ''' import numpy as np ds = [360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532] ds = np.array(ds)/1e6 numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1] dist = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0) psd = PSDInterpolated(dist.Dis, dist.fractions) assert len(psd.fractions) == len(psd.ds) assert len(psd.fractions) == 15 # test fractions_discrete vs input assert_close1d(psd.fractions_discrete(ds), dist.fractions) # test cdf_discrete assert_close1d(psd.cdf_discrete(ds), psd.fraction_cdf[1:]) # test that dn solves backwards for exactly the right value - slow cumulative_fractions = np.cumsum(dist.fractions) ds_for_fractions = np.array([psd.dn(f) for f in cumulative_fractions]) assert_close1d(ds, ds_for_fractions) # test _pdf test_pdf = psd._pdf(1e-3) assert_close(test_pdf, 106.28284463095554) # test _cdf test_cdf = psd._cdf(1e-3) assert_close(test_cdf, 0.02278897476363087) # test _pdf_basis_integral test_int = psd._pdf_basis_integral(1e-3, 2) assert_close(test_int, 1.509707233427664e-08) # Check that the 0 point is created and the points and fractions are the same assert_close1d(psd.ds, [0] + ds.tolist()) assert_close1d(psd.fractions, [0] + dist.fractions) # test mean_size test_mean = psd.mean_size(3, 2) assert_close(test_mean, 0.002211577679574544) def test_PSDInterpolated_discrete_range_pts(): import numpy as np ds = [360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532] ds = [i*1e-6 for i in ds] numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1] psd = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0) # test fractions_discrete vs input assert_close1d(psd.fractions_discrete(ds), psd.fractions) # test cdf_discrete assert_close1d(psd.cdf_discrete(ds), psd.interpolated.fraction_cdf[1:]) # test that dn solves backwards for exactly the right value cumulative_fractions = np.cumsum(psd.fractions) ds_for_fractions = np.array([psd.dn(f) for f in cumulative_fractions]) assert_close1d(ds, ds_for_fractions) # test _pdf test_pdf = psd.pdf(1e-3) assert_close(test_pdf, 106.28284463095554) # test _cdf test_cdf = psd.cdf(1e-3) assert_close(test_cdf, 0.02278897476363087) # test _pdf_basis_integral test_int = psd._pdf_basis_integral(1e-3, 2) assert_close(test_int, 1.509707233427664e-08) assert not isclose(psd.mean_size(3, 2), psd.interpolated.mean_size(3, 2)) def test_PSDInterpolated_discrete(): import numpy as np ds = 1E-6*np.array([240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532]) numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1] psd = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0) dn_05 = psd.dn(0.5) assert_close(dn_05, 0.002526452452632658) assert_close(psd.cdf(dn_05), 0.5, rtol=1e-4) assert_close(psd.cdf(0.002, 4), 0.18010145594167873, rtol=5e-3) assert_close(psd.pdf(0.002, 2), 466.42761174380007, rtol=5e-3) assert_close(psd.pdf(0.002526452452632658), 408.4774528377241, rtol=1e-2) assert_close(psd.mean_size(1, 0), psd.interpolated.mean_size(1, 0), rtol=0.05) assert_close(psd.cdf(100), 1) assert_close(psd.cdf(psd.dn(1)), 1) def test_psd_spacing(): ans_log = psd_spacing(d_min=1, d_max=10, pts=4, method='logarithmic') ans_log_expect = [1.0, 2.154434690031884, 4.641588833612778, 10.0] assert_close1d(ans_log, ans_log_expect) ans_lin = psd_spacing(d_min=0, d_max=10, pts=4, method='linear') ans_lin_expect = [0.0, 3.3333333333333335, 6.666666666666667, 10.0] assert_close1d(ans_lin, ans_lin_expect) with pytest.raises(Exception): psd_spacing(d_min=0, d_max=10, pts=8, method='R5') with pytest.raises(Exception): psd_spacing(d_min=5e-5, d_max=5e-4, method='BADMETHOD') # This example from an iso standard, ISO 9276-2 2014 ans_R5 = psd_spacing(d_max=25, pts=8, method='R5') ans_R5_expect = [0.9952679263837426, 1.5773933612004825, 2.499999999999999, 3.9622329811527823, 6.279716078773949, 9.95267926383743, 15.77393361200483, 25] assert_close1d(ans_R5, ans_R5_expect) ans_R5_reversed = psd_spacing(d_min=0.9952679263837426, pts=8, method='R5') assert_close1d(ans_R5_reversed, ans_R5_expect) ans_R5_float = psd_spacing(d_max=25, pts=8, method='R5.00000001') assert_close1d(ans_R5_float, ans_R5_expect) ans = psd_spacing(d_min=5e-5, d_max=1e-3, method='ISO 3310-1') ans_expect = [5e-05, 5.3e-05, 5.6e-05, 6.3e-05, 7.1e-05, 7.5e-05, 8e-05, 9e-05, 0.0001, 0.000106, 0.000112, 0.000125, 0.00014, 0.00015, 0.00016, 0.00018, 0.0002, 0.000212, 0.000224, 0.00025, 0.00028, 0.0003, 0.000315, 0.000355, 0.0004, 0.000425, 0.00045, 0.0005, 0.00056, 0.0006, 0.00063, 0.00071, 0.0008, 0.00085, 0.0009, 0.001] assert_close1d(ans, ans_expect) assert [] == psd_spacing(d_min=0, d_max=1e-6, method='ISO 3310-1') assert [] == psd_spacing(d_min=1, d_max=1e2, method='ISO 3310-1') assert psd_spacing(d_min=5e-5, d_max=1e-3, method='ISO 3310-1 R20') assert psd_spacing(d_min=5e-5, d_max=1e-3, method='ISO 3310-1 R20/3') assert psd_spacing(d_min=5e-5, d_max=1e-3, method='ISO 3310-1 R40/3') assert psd_spacing(d_min=0e-5, d_max=1e-3, method='ISO 3310-1 R10') ds = psd_spacing(d_min=1e-5, d_max=1e-4, method='ASTM E11') ds_expect = [2e-05, 2.5e-05, 3.2e-05, 3.8e-05, 4.5e-05, 5.3e-05, 6.3e-05, 7.5e-05, 9e-05] assert_close1d(ds, ds_expect) def test_example_integer_power_issue(): diameters = [2, 5, 10, 15, 20, 30, 40, 50] # Particle diameters in micrometers (µm) volume_percentages = [5, 10, 20, 25, 20, 10, 5, 5] # Volume % in each size class cement_psd = ParticleSizeDistribution(ds=diameters, fractions=volume_percentages, order=3) val = cement_psd.pdf(4.61538459884405e-08, n=2) assert val is not None def test_PSD_numpy_input(): import numpy as np ds = np.array([240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532]) numbers = np.array([65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1]) psd = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0) assert type(psd.mean_size(2, 1)) is float assert type(psd.mean_size(1, 0)) is float assert type(psd.mean_size(3, 2)) is float assert type(psd.mean_size(4, 3)) is float assert type(psd.dn(0.1)) is float assert type(psd.dn(0.9)) is float assert type(psd.pdf(1000)) is float assert type(psd.cdf(5000)) is float psd = PSDLognormal(s=0.5, d_characteristic=5E-6) assert type(psd.pdf(1e-6)) is float assert type(psd.cdf(7e-6)) is float assert type(psd.dn(0.1)) is float assert type(psd.mean_size(3, 2)) is float ds_discrete = psd.ds_discrete(pts=8) fractions_discrete = psd.fractions_discrete(ds_discrete) psd_discrete = ParticleSizeDistribution(ds=ds_discrete, fractions=fractions_discrete, order=3) assert type(psd_discrete.mean_size(3, 2)) is float assert type(psd.vssa) is float assert type(psd.dn(.9) - psd.dn(0.1)) is float assert type((psd.dn(.9) - psd.dn(0.1))/psd.dn(0.5)) is float assert type(psd.dn(0.75)/psd.dn(0.25)) is float assert type(psd.dn(0.9)/psd.dn(0.1)) is float