1# Created by Pearu Peterson, June 2003 2import numpy as np 3from numpy.testing import (assert_equal, assert_almost_equal, assert_array_equal, 4 assert_array_almost_equal, assert_allclose, suppress_warnings) 5from pytest import raises as assert_raises 6 7from numpy import array, diff, linspace, meshgrid, ones, pi, shape 8from scipy.interpolate.fitpack import bisplrep, bisplev 9from scipy.interpolate.fitpack2 import (UnivariateSpline, 10 LSQUnivariateSpline, InterpolatedUnivariateSpline, 11 LSQBivariateSpline, SmoothBivariateSpline, RectBivariateSpline, 12 LSQSphereBivariateSpline, SmoothSphereBivariateSpline, 13 RectSphereBivariateSpline) 14 15 16class TestUnivariateSpline: 17 def test_linear_constant(self): 18 x = [1,2,3] 19 y = [3,3,3] 20 lut = UnivariateSpline(x,y,k=1) 21 assert_array_almost_equal(lut.get_knots(),[1,3]) 22 assert_array_almost_equal(lut.get_coeffs(),[3,3]) 23 assert_almost_equal(lut.get_residual(),0.0) 24 assert_array_almost_equal(lut([1,1.5,2]),[3,3,3]) 25 26 def test_preserve_shape(self): 27 x = [1, 2, 3] 28 y = [0, 2, 4] 29 lut = UnivariateSpline(x, y, k=1) 30 arg = 2 31 assert_equal(shape(arg), shape(lut(arg))) 32 assert_equal(shape(arg), shape(lut(arg, nu=1))) 33 arg = [1.5, 2, 2.5] 34 assert_equal(shape(arg), shape(lut(arg))) 35 assert_equal(shape(arg), shape(lut(arg, nu=1))) 36 37 def test_linear_1d(self): 38 x = [1,2,3] 39 y = [0,2,4] 40 lut = UnivariateSpline(x,y,k=1) 41 assert_array_almost_equal(lut.get_knots(),[1,3]) 42 assert_array_almost_equal(lut.get_coeffs(),[0,4]) 43 assert_almost_equal(lut.get_residual(),0.0) 44 assert_array_almost_equal(lut([1,1.5,2]),[0,1,2]) 45 46 def test_subclassing(self): 47 # See #731 48 49 class ZeroSpline(UnivariateSpline): 50 def __call__(self, x): 51 return 0*array(x) 52 53 sp = ZeroSpline([1,2,3,4,5], [3,2,3,2,3], k=2) 54 assert_array_equal(sp([1.5, 2.5]), [0., 0.]) 55 56 def test_empty_input(self): 57 # Test whether empty input returns an empty output. Ticket 1014 58 x = [1,3,5,7,9] 59 y = [0,4,9,12,21] 60 spl = UnivariateSpline(x, y, k=3) 61 assert_array_equal(spl([]), array([])) 62 63 def test_resize_regression(self): 64 """Regression test for #1375.""" 65 x = [-1., -0.65016502, -0.58856235, -0.26903553, -0.17370892, 66 -0.10011001, 0., 0.10011001, 0.17370892, 0.26903553, 0.58856235, 67 0.65016502, 1.] 68 y = [1.,0.62928599, 0.5797223, 0.39965815, 0.36322694, 0.3508061, 69 0.35214793, 0.3508061, 0.36322694, 0.39965815, 0.5797223, 70 0.62928599, 1.] 71 w = [1.00000000e+12, 6.88875973e+02, 4.89314737e+02, 4.26864807e+02, 72 6.07746770e+02, 4.51341444e+02, 3.17480210e+02, 4.51341444e+02, 73 6.07746770e+02, 4.26864807e+02, 4.89314737e+02, 6.88875973e+02, 74 1.00000000e+12] 75 spl = UnivariateSpline(x=x, y=y, w=w, s=None) 76 desired = array([0.35100374, 0.51715855, 0.87789547, 0.98719344]) 77 assert_allclose(spl([0.1, 0.5, 0.9, 0.99]), desired, atol=5e-4) 78 79 def test_out_of_range_regression(self): 80 # Test different extrapolation modes. See ticket 3557 81 x = np.arange(5, dtype=float) 82 y = x**3 83 84 xp = linspace(-8, 13, 100) 85 xp_zeros = xp.copy() 86 xp_zeros[np.logical_or(xp_zeros < 0., xp_zeros > 4.)] = 0 87 xp_clip = xp.copy() 88 xp_clip[xp_clip < x[0]] = x[0] 89 xp_clip[xp_clip > x[-1]] = x[-1] 90 91 for cls in [UnivariateSpline, InterpolatedUnivariateSpline]: 92 spl = cls(x=x, y=y) 93 for ext in [0, 'extrapolate']: 94 assert_allclose(spl(xp, ext=ext), xp**3, atol=1e-16) 95 assert_allclose(cls(x, y, ext=ext)(xp), xp**3, atol=1e-16) 96 for ext in [1, 'zeros']: 97 assert_allclose(spl(xp, ext=ext), xp_zeros**3, atol=1e-16) 98 assert_allclose(cls(x, y, ext=ext)(xp), xp_zeros**3, atol=1e-16) 99 for ext in [2, 'raise']: 100 assert_raises(ValueError, spl, xp, **dict(ext=ext)) 101 for ext in [3, 'const']: 102 assert_allclose(spl(xp, ext=ext), xp_clip**3, atol=1e-16) 103 assert_allclose(cls(x, y, ext=ext)(xp), xp_clip**3, atol=1e-16) 104 105 # also test LSQUnivariateSpline [which needs explicit knots] 106 t = spl.get_knots()[3:4] # interior knots w/ default k=3 107 spl = LSQUnivariateSpline(x, y, t) 108 assert_allclose(spl(xp, ext=0), xp**3, atol=1e-16) 109 assert_allclose(spl(xp, ext=1), xp_zeros**3, atol=1e-16) 110 assert_raises(ValueError, spl, xp, **dict(ext=2)) 111 assert_allclose(spl(xp, ext=3), xp_clip**3, atol=1e-16) 112 113 # also make sure that unknown values for `ext` are caught early 114 for ext in [-1, 'unknown']: 115 spl = UnivariateSpline(x, y) 116 assert_raises(ValueError, spl, xp, **dict(ext=ext)) 117 assert_raises(ValueError, UnivariateSpline, 118 **dict(x=x, y=y, ext=ext)) 119 120 def test_lsq_fpchec(self): 121 xs = np.arange(100) * 1. 122 ys = np.arange(100) * 1. 123 knots = np.linspace(0, 99, 10) 124 bbox = (-1, 101) 125 assert_raises(ValueError, LSQUnivariateSpline, xs, ys, knots, 126 bbox=bbox) 127 128 def test_derivative_and_antiderivative(self): 129 # Thin wrappers to splder/splantider, so light smoke test only. 130 x = np.linspace(0, 1, 70)**3 131 y = np.cos(x) 132 133 spl = UnivariateSpline(x, y, s=0) 134 spl2 = spl.antiderivative(2).derivative(2) 135 assert_allclose(spl(0.3), spl2(0.3)) 136 137 spl2 = spl.antiderivative(1) 138 assert_allclose(spl2(0.6) - spl2(0.2), 139 spl.integral(0.2, 0.6)) 140 141 def test_derivative_extrapolation(self): 142 # Regression test for gh-10195: for a const-extrapolation spline 143 # its derivative evaluates to zero for extrapolation 144 x_values = [1, 2, 4, 6, 8.5] 145 y_values = [0.5, 0.8, 1.3, 2.5, 5] 146 f = UnivariateSpline(x_values, y_values, ext='const', k=3) 147 148 x = [-1, 0, -0.5, 9, 9.5, 10] 149 assert_allclose(f.derivative()(x), 0, atol=1e-15) 150 151 def test_integral_out_of_bounds(self): 152 # Regression test for gh-7906: .integral(a, b) is wrong if both 153 # a and b are out-of-bounds 154 x = np.linspace(0., 1., 7) 155 for ext in range(4): 156 f = UnivariateSpline(x, x, s=0, ext=ext) 157 for (a, b) in [(1, 1), (1, 5), (2, 5), 158 (0, 0), (-2, 0), (-2, -1)]: 159 assert_allclose(f.integral(a, b), 0, atol=1e-15) 160 161 def test_nan(self): 162 # bail out early if the input data contains nans 163 x = np.arange(10, dtype=float) 164 y = x**3 165 w = np.ones_like(x) 166 # also test LSQUnivariateSpline [which needs explicit knots] 167 spl = UnivariateSpline(x, y, check_finite=True) 168 t = spl.get_knots()[3:4] # interior knots w/ default k=3 169 y_end = y[-1] 170 for z in [np.nan, np.inf, -np.inf]: 171 y[-1] = z 172 assert_raises(ValueError, UnivariateSpline, 173 **dict(x=x, y=y, check_finite=True)) 174 assert_raises(ValueError, InterpolatedUnivariateSpline, 175 **dict(x=x, y=y, check_finite=True)) 176 assert_raises(ValueError, LSQUnivariateSpline, 177 **dict(x=x, y=y, t=t, check_finite=True)) 178 y[-1] = y_end # check valid y but invalid w 179 w[-1] = z 180 assert_raises(ValueError, UnivariateSpline, 181 **dict(x=x, y=y, w=w, check_finite=True)) 182 assert_raises(ValueError, InterpolatedUnivariateSpline, 183 **dict(x=x, y=y, w=w, check_finite=True)) 184 assert_raises(ValueError, LSQUnivariateSpline, 185 **dict(x=x, y=y, t=t, w=w, check_finite=True)) 186 187 def test_strictly_increasing_x(self): 188 # Test the x is required to be strictly increasing for 189 # UnivariateSpline if s=0 and for InterpolatedUnivariateSpline, 190 # but merely increasing for UnivariateSpline if s>0 191 # and for LSQUnivariateSpline; see gh-8535 192 xx = np.arange(10, dtype=float) 193 yy = xx**3 194 x = np.arange(10, dtype=float) 195 x[1] = x[0] 196 y = x**3 197 w = np.ones_like(x) 198 # also test LSQUnivariateSpline [which needs explicit knots] 199 spl = UnivariateSpline(xx, yy, check_finite=True) 200 t = spl.get_knots()[3:4] # interior knots w/ default k=3 201 UnivariateSpline(x=x, y=y, w=w, s=1, check_finite=True) 202 LSQUnivariateSpline(x=x, y=y, t=t, w=w, check_finite=True) 203 assert_raises(ValueError, UnivariateSpline, 204 **dict(x=x, y=y, s=0, check_finite=True)) 205 assert_raises(ValueError, InterpolatedUnivariateSpline, 206 **dict(x=x, y=y, check_finite=True)) 207 208 def test_increasing_x(self): 209 # Test that x is required to be increasing, see gh-8535 210 xx = np.arange(10, dtype=float) 211 yy = xx**3 212 x = np.arange(10, dtype=float) 213 x[1] = x[0] - 1.0 214 y = x**3 215 w = np.ones_like(x) 216 # also test LSQUnivariateSpline [which needs explicit knots] 217 spl = UnivariateSpline(xx, yy, check_finite=True) 218 t = spl.get_knots()[3:4] # interior knots w/ default k=3 219 assert_raises(ValueError, UnivariateSpline, 220 **dict(x=x, y=y, check_finite=True)) 221 assert_raises(ValueError, InterpolatedUnivariateSpline, 222 **dict(x=x, y=y, check_finite=True)) 223 assert_raises(ValueError, LSQUnivariateSpline, 224 **dict(x=x, y=y, t=t, w=w, check_finite=True)) 225 226 def test_invalid_input_for_univariate_spline(self): 227 228 with assert_raises(ValueError) as info: 229 x_values = [1, 2, 4, 6, 8.5] 230 y_values = [0.5, 0.8, 1.3, 2.5] 231 UnivariateSpline(x_values, y_values) 232 assert "x and y should have a same length" in str(info.value) 233 234 with assert_raises(ValueError) as info: 235 x_values = [1, 2, 4, 6, 8.5] 236 y_values = [0.5, 0.8, 1.3, 2.5, 2.8] 237 w_values = [-1.0, 1.0, 1.0, 1.0] 238 UnivariateSpline(x_values, y_values, w=w_values) 239 assert "x, y, and w should have a same length" in str(info.value) 240 241 with assert_raises(ValueError) as info: 242 bbox = (-1) 243 UnivariateSpline(x_values, y_values, bbox=bbox) 244 assert "bbox shape should be (2,)" in str(info.value) 245 246 with assert_raises(ValueError) as info: 247 UnivariateSpline(x_values, y_values, k=6) 248 assert "k should be 1 <= k <= 5" in str(info.value) 249 250 with assert_raises(ValueError) as info: 251 UnivariateSpline(x_values, y_values, s=-1.0) 252 assert "s should be s >= 0.0" in str(info.value) 253 254 def test_invalid_input_for_interpolated_univariate_spline(self): 255 256 with assert_raises(ValueError) as info: 257 x_values = [1, 2, 4, 6, 8.5] 258 y_values = [0.5, 0.8, 1.3, 2.5] 259 InterpolatedUnivariateSpline(x_values, y_values) 260 assert "x and y should have a same length" in str(info.value) 261 262 with assert_raises(ValueError) as info: 263 x_values = [1, 2, 4, 6, 8.5] 264 y_values = [0.5, 0.8, 1.3, 2.5, 2.8] 265 w_values = [-1.0, 1.0, 1.0, 1.0] 266 InterpolatedUnivariateSpline(x_values, y_values, w=w_values) 267 assert "x, y, and w should have a same length" in str(info.value) 268 269 with assert_raises(ValueError) as info: 270 bbox = (-1) 271 InterpolatedUnivariateSpline(x_values, y_values, bbox=bbox) 272 assert "bbox shape should be (2,)" in str(info.value) 273 274 with assert_raises(ValueError) as info: 275 InterpolatedUnivariateSpline(x_values, y_values, k=6) 276 assert "k should be 1 <= k <= 5" in str(info.value) 277 278 def test_invalid_input_for_lsq_univariate_spline(self): 279 280 x_values = [1, 2, 4, 6, 8.5] 281 y_values = [0.5, 0.8, 1.3, 2.5, 2.8] 282 spl = UnivariateSpline(x_values, y_values, check_finite=True) 283 t_values = spl.get_knots()[3:4] # interior knots w/ default k=3 284 285 with assert_raises(ValueError) as info: 286 x_values = [1, 2, 4, 6, 8.5] 287 y_values = [0.5, 0.8, 1.3, 2.5] 288 LSQUnivariateSpline(x_values, y_values, t_values) 289 assert "x and y should have a same length" in str(info.value) 290 291 with assert_raises(ValueError) as info: 292 x_values = [1, 2, 4, 6, 8.5] 293 y_values = [0.5, 0.8, 1.3, 2.5, 2.8] 294 w_values = [1.0, 1.0, 1.0, 1.0] 295 LSQUnivariateSpline(x_values, y_values, t_values, w=w_values) 296 assert "x, y, and w should have a same length" in str(info.value) 297 298 with assert_raises(ValueError) as info: 299 bbox = (100, -100) 300 LSQUnivariateSpline(x_values, y_values, t_values, bbox=bbox) 301 assert "Interior knots t must satisfy Schoenberg-Whitney conditions" in str(info.value) 302 303 with assert_raises(ValueError) as info: 304 bbox = (-1) 305 LSQUnivariateSpline(x_values, y_values, t_values, bbox=bbox) 306 assert "bbox shape should be (2,)" in str(info.value) 307 308 with assert_raises(ValueError) as info: 309 LSQUnivariateSpline(x_values, y_values, t_values, k=6) 310 assert "k should be 1 <= k <= 5" in str(info.value) 311 312 def test_array_like_input(self): 313 x_values = np.array([1, 2, 4, 6, 8.5]) 314 y_values = np.array([0.5, 0.8, 1.3, 2.5, 2.8]) 315 w_values = np.array([1.0, 1.0, 1.0, 1.0, 1.0]) 316 bbox = np.array([-100, 100]) 317 # np.array input 318 spl1 = UnivariateSpline(x=x_values, y=y_values, w=w_values, 319 bbox=bbox) 320 # list input 321 spl2 = UnivariateSpline(x=x_values.tolist(), y=y_values.tolist(), 322 w=w_values.tolist(), bbox=bbox.tolist()) 323 324 assert_allclose(spl1([0.1, 0.5, 0.9, 0.99]), 325 spl2([0.1, 0.5, 0.9, 0.99])) 326 327 328class TestLSQBivariateSpline: 329 # NOTE: The systems in this test class are rank-deficient 330 def test_linear_constant(self): 331 x = [1,1,1,2,2,2,3,3,3] 332 y = [1,2,3,1,2,3,1,2,3] 333 z = [3,3,3,3,3,3,3,3,3] 334 s = 0.1 335 tx = [1+s,3-s] 336 ty = [1+s,3-s] 337 with suppress_warnings() as sup: 338 r = sup.record(UserWarning, "\nThe coefficients of the spline") 339 lut = LSQBivariateSpline(x,y,z,tx,ty,kx=1,ky=1) 340 assert_equal(len(r), 1) 341 342 assert_almost_equal(lut(2,2), 3.) 343 344 def test_bilinearity(self): 345 x = [1,1,1,2,2,2,3,3,3] 346 y = [1,2,3,1,2,3,1,2,3] 347 z = [0,7,8,3,4,7,1,3,4] 348 s = 0.1 349 tx = [1+s,3-s] 350 ty = [1+s,3-s] 351 with suppress_warnings() as sup: 352 # This seems to fail (ier=1, see ticket 1642). 353 sup.filter(UserWarning, "\nThe coefficients of the spline") 354 lut = LSQBivariateSpline(x,y,z,tx,ty,kx=1,ky=1) 355 356 tx, ty = lut.get_knots() 357 for xa, xb in zip(tx[:-1], tx[1:]): 358 for ya, yb in zip(ty[:-1], ty[1:]): 359 for t in [0.1, 0.5, 0.9]: 360 for s in [0.3, 0.4, 0.7]: 361 xp = xa*(1-t) + xb*t 362 yp = ya*(1-s) + yb*s 363 zp = (+ lut(xa, ya)*(1-t)*(1-s) 364 + lut(xb, ya)*t*(1-s) 365 + lut(xa, yb)*(1-t)*s 366 + lut(xb, yb)*t*s) 367 assert_almost_equal(lut(xp,yp), zp) 368 369 def test_integral(self): 370 x = [1,1,1,2,2,2,8,8,8] 371 y = [1,2,3,1,2,3,1,2,3] 372 z = array([0,7,8,3,4,7,1,3,4]) 373 374 s = 0.1 375 tx = [1+s,3-s] 376 ty = [1+s,3-s] 377 with suppress_warnings() as sup: 378 r = sup.record(UserWarning, "\nThe coefficients of the spline") 379 lut = LSQBivariateSpline(x, y, z, tx, ty, kx=1, ky=1) 380 assert_equal(len(r), 1) 381 tx, ty = lut.get_knots() 382 tz = lut(tx, ty) 383 trpz = .25*(diff(tx)[:,None]*diff(ty)[None,:] 384 * (tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum() 385 386 assert_almost_equal(lut.integral(tx[0], tx[-1], ty[0], ty[-1]), 387 trpz) 388 389 def test_empty_input(self): 390 # Test whether empty inputs returns an empty output. Ticket 1014 391 x = [1,1,1,2,2,2,3,3,3] 392 y = [1,2,3,1,2,3,1,2,3] 393 z = [3,3,3,3,3,3,3,3,3] 394 s = 0.1 395 tx = [1+s,3-s] 396 ty = [1+s,3-s] 397 with suppress_warnings() as sup: 398 r = sup.record(UserWarning, "\nThe coefficients of the spline") 399 lut = LSQBivariateSpline(x, y, z, tx, ty, kx=1, ky=1) 400 assert_equal(len(r), 1) 401 402 assert_array_equal(lut([], []), np.zeros((0,0))) 403 assert_array_equal(lut([], [], grid=False), np.zeros((0,))) 404 405 def test_invalid_input(self): 406 s = 0.1 407 tx = [1 + s, 3 - s] 408 ty = [1 + s, 3 - s] 409 410 with assert_raises(ValueError) as info: 411 x = np.linspace(1.0, 10.0) 412 y = np.linspace(1.0, 10.0) 413 z = np.linspace(1.0, 10.0, num=10) 414 LSQBivariateSpline(x, y, z, tx, ty) 415 assert "x, y, and z should have a same length" in str(info.value) 416 417 with assert_raises(ValueError) as info: 418 x = np.linspace(1.0, 10.0) 419 y = np.linspace(1.0, 10.0) 420 z = np.linspace(1.0, 10.0) 421 w = np.linspace(1.0, 10.0, num=20) 422 LSQBivariateSpline(x, y, z, tx, ty, w=w) 423 assert "x, y, z, and w should have a same length" in str(info.value) 424 425 with assert_raises(ValueError) as info: 426 w = np.linspace(-1.0, 10.0) 427 LSQBivariateSpline(x, y, z, tx, ty, w=w) 428 assert "w should be positive" in str(info.value) 429 430 with assert_raises(ValueError) as info: 431 bbox = (-100, 100, -100) 432 LSQBivariateSpline(x, y, z, tx, ty, bbox=bbox) 433 assert "bbox shape should be (4,)" in str(info.value) 434 435 with assert_raises(ValueError) as info: 436 LSQBivariateSpline(x, y, z, tx, ty, kx=10, ky=10) 437 assert "The length of x, y and z should be at least (kx+1) * (ky+1)" in \ 438 str(info.value) 439 440 with assert_raises(ValueError) as exc_info: 441 LSQBivariateSpline(x, y, z, tx, ty, eps=0.0) 442 assert "eps should be between (0, 1)" in str(exc_info.value) 443 444 with assert_raises(ValueError) as exc_info: 445 LSQBivariateSpline(x, y, z, tx, ty, eps=1.0) 446 assert "eps should be between (0, 1)" in str(exc_info.value) 447 448 def test_array_like_input(self): 449 s = 0.1 450 tx = np.array([1 + s, 3 - s]) 451 ty = np.array([1 + s, 3 - s]) 452 x = np.linspace(1.0, 10.0) 453 y = np.linspace(1.0, 10.0) 454 z = np.linspace(1.0, 10.0) 455 w = np.linspace(1.0, 10.0) 456 bbox = np.array([1.0, 10.0, 1.0, 10.0]) 457 458 with suppress_warnings() as sup: 459 r = sup.record(UserWarning, "\nThe coefficients of the spline") 460 # np.array input 461 spl1 = LSQBivariateSpline(x, y, z, tx, ty, w=w, bbox=bbox) 462 # list input 463 spl2 = LSQBivariateSpline(x.tolist(), y.tolist(), z.tolist(), 464 tx.tolist(), ty.tolist(), w=w.tolist(), 465 bbox=bbox) 466 assert_allclose(spl1(2.0, 2.0), spl2(2.0, 2.0)) 467 assert_equal(len(r), 2) 468 469 def test_unequal_length_of_knots(self): 470 """Test for the case when the input knot-location arrays in x and y are 471 of different lengths. 472 """ 473 x, y = np.mgrid[0:100, 0:100] 474 x = x.ravel() 475 y = y.ravel() 476 z = 3.0 * np.ones_like(x) 477 tx = np.linspace(0.1, 98.0, 29) 478 ty = np.linspace(0.1, 98.0, 33) 479 with suppress_warnings() as sup: 480 r = sup.record(UserWarning, "\nThe coefficients of the spline") 481 lut = LSQBivariateSpline(x,y,z,tx,ty) 482 assert_equal(len(r), 1) 483 484 assert_almost_equal(lut(x, y, grid=False), z) 485 486 487class TestSmoothBivariateSpline: 488 def test_linear_constant(self): 489 x = [1,1,1,2,2,2,3,3,3] 490 y = [1,2,3,1,2,3,1,2,3] 491 z = [3,3,3,3,3,3,3,3,3] 492 lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1) 493 assert_array_almost_equal(lut.get_knots(),([1,1,3,3],[1,1,3,3])) 494 assert_array_almost_equal(lut.get_coeffs(),[3,3,3,3]) 495 assert_almost_equal(lut.get_residual(),0.0) 496 assert_array_almost_equal(lut([1,1.5,2],[1,1.5]),[[3,3],[3,3],[3,3]]) 497 498 def test_linear_1d(self): 499 x = [1,1,1,2,2,2,3,3,3] 500 y = [1,2,3,1,2,3,1,2,3] 501 z = [0,0,0,2,2,2,4,4,4] 502 lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1) 503 assert_array_almost_equal(lut.get_knots(),([1,1,3,3],[1,1,3,3])) 504 assert_array_almost_equal(lut.get_coeffs(),[0,0,4,4]) 505 assert_almost_equal(lut.get_residual(),0.0) 506 assert_array_almost_equal(lut([1,1.5,2],[1,1.5]),[[0,0],[1,1],[2,2]]) 507 508 def test_integral(self): 509 x = [1,1,1,2,2,2,4,4,4] 510 y = [1,2,3,1,2,3,1,2,3] 511 z = array([0,7,8,3,4,7,1,3,4]) 512 513 with suppress_warnings() as sup: 514 # This seems to fail (ier=1, see ticket 1642). 515 sup.filter(UserWarning, "\nThe required storage space") 516 lut = SmoothBivariateSpline(x, y, z, kx=1, ky=1, s=0) 517 518 tx = [1,2,4] 519 ty = [1,2,3] 520 521 tz = lut(tx, ty) 522 trpz = .25*(diff(tx)[:,None]*diff(ty)[None,:] 523 * (tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum() 524 assert_almost_equal(lut.integral(tx[0], tx[-1], ty[0], ty[-1]), trpz) 525 526 lut2 = SmoothBivariateSpline(x, y, z, kx=2, ky=2, s=0) 527 assert_almost_equal(lut2.integral(tx[0], tx[-1], ty[0], ty[-1]), trpz, 528 decimal=0) # the quadratures give 23.75 and 23.85 529 530 tz = lut(tx[:-1], ty[:-1]) 531 trpz = .25*(diff(tx[:-1])[:,None]*diff(ty[:-1])[None,:] 532 * (tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum() 533 assert_almost_equal(lut.integral(tx[0], tx[-2], ty[0], ty[-2]), trpz) 534 535 def test_rerun_lwrk2_too_small(self): 536 # in this setting, lwrk2 is too small in the default run. Here we 537 # check for equality with the bisplrep/bisplev output because there, 538 # an automatic re-run of the spline representation is done if ier>10. 539 x = np.linspace(-2, 2, 80) 540 y = np.linspace(-2, 2, 80) 541 z = x + y 542 xi = np.linspace(-1, 1, 100) 543 yi = np.linspace(-2, 2, 100) 544 tck = bisplrep(x, y, z) 545 res1 = bisplev(xi, yi, tck) 546 interp_ = SmoothBivariateSpline(x, y, z) 547 res2 = interp_(xi, yi) 548 assert_almost_equal(res1, res2) 549 550 def test_invalid_input(self): 551 552 with assert_raises(ValueError) as info: 553 x = np.linspace(1.0, 10.0) 554 y = np.linspace(1.0, 10.0) 555 z = np.linspace(1.0, 10.0, num=10) 556 SmoothBivariateSpline(x, y, z) 557 assert "x, y, and z should have a same length" in str(info.value) 558 559 with assert_raises(ValueError) as info: 560 x = np.linspace(1.0, 10.0) 561 y = np.linspace(1.0, 10.0) 562 z = np.linspace(1.0, 10.0) 563 w = np.linspace(1.0, 10.0, num=20) 564 SmoothBivariateSpline(x, y, z, w=w) 565 assert "x, y, z, and w should have a same length" in str(info.value) 566 567 with assert_raises(ValueError) as info: 568 w = np.linspace(-1.0, 10.0) 569 SmoothBivariateSpline(x, y, z, w=w) 570 assert "w should be positive" in str(info.value) 571 572 with assert_raises(ValueError) as info: 573 bbox = (-100, 100, -100) 574 SmoothBivariateSpline(x, y, z, bbox=bbox) 575 assert "bbox shape should be (4,)" in str(info.value) 576 577 with assert_raises(ValueError) as info: 578 SmoothBivariateSpline(x, y, z, kx=10, ky=10) 579 assert "The length of x, y and z should be at least (kx+1) * (ky+1)" in\ 580 str(info.value) 581 582 with assert_raises(ValueError) as info: 583 SmoothBivariateSpline(x, y, z, s=-1.0) 584 assert "s should be s >= 0.0" in str(info.value) 585 586 with assert_raises(ValueError) as exc_info: 587 SmoothBivariateSpline(x, y, z, eps=0.0) 588 assert "eps should be between (0, 1)" in str(exc_info.value) 589 590 with assert_raises(ValueError) as exc_info: 591 SmoothBivariateSpline(x, y, z, eps=1.0) 592 assert "eps should be between (0, 1)" in str(exc_info.value) 593 594 def test_array_like_input(self): 595 x = np.array([1, 1, 1, 2, 2, 2, 3, 3, 3]) 596 y = np.array([1, 2, 3, 1, 2, 3, 1, 2, 3]) 597 z = np.array([3, 3, 3, 3, 3, 3, 3, 3, 3]) 598 w = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1]) 599 bbox = np.array([1.0, 3.0, 1.0, 3.0]) 600 # np.array input 601 spl1 = SmoothBivariateSpline(x, y, z, w=w, bbox=bbox, kx=1, ky=1) 602 # list input 603 spl2 = SmoothBivariateSpline(x.tolist(), y.tolist(), z.tolist(), 604 bbox=bbox.tolist(), w=w.tolist(), 605 kx=1, ky=1) 606 assert_allclose(spl1(0.1, 0.5), spl2(0.1, 0.5)) 607 608 609class TestLSQSphereBivariateSpline: 610 def setup_method(self): 611 # define the input data and coordinates 612 ntheta, nphi = 70, 90 613 theta = linspace(0.5/(ntheta - 1), 1 - 0.5/(ntheta - 1), ntheta) * pi 614 phi = linspace(0.5/(nphi - 1), 1 - 0.5/(nphi - 1), nphi) * 2. * pi 615 data = ones((theta.shape[0], phi.shape[0])) 616 # define knots and extract data values at the knots 617 knotst = theta[::5] 618 knotsp = phi[::5] 619 knotdata = data[::5, ::5] 620 # calculate spline coefficients 621 lats, lons = meshgrid(theta, phi) 622 lut_lsq = LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), 623 data.T.ravel(), knotst, knotsp) 624 self.lut_lsq = lut_lsq 625 self.data = knotdata 626 self.new_lons, self.new_lats = knotsp, knotst 627 628 def test_linear_constant(self): 629 assert_almost_equal(self.lut_lsq.get_residual(), 0.0) 630 assert_array_almost_equal(self.lut_lsq(self.new_lats, self.new_lons), 631 self.data) 632 633 def test_empty_input(self): 634 assert_array_almost_equal(self.lut_lsq([], []), np.zeros((0,0))) 635 assert_array_almost_equal(self.lut_lsq([], [], grid=False), np.zeros((0,))) 636 637 def test_invalid_input(self): 638 ntheta, nphi = 70, 90 639 theta = linspace(0.5 / (ntheta - 1), 1 - 0.5 / (ntheta - 1), 640 ntheta) * pi 641 phi = linspace(0.5 / (nphi - 1), 1 - 0.5 / (nphi - 1), nphi) * 2. * pi 642 data = ones((theta.shape[0], phi.shape[0])) 643 # define knots and extract data values at the knots 644 knotst = theta[::5] 645 knotsp = phi[::5] 646 647 with assert_raises(ValueError) as exc_info: 648 invalid_theta = linspace(-0.1, 1.0, num=ntheta) * pi 649 invalid_lats, lons = meshgrid(invalid_theta, phi) 650 LSQSphereBivariateSpline(invalid_lats.ravel(), lons.ravel(), 651 data.T.ravel(), knotst, knotsp) 652 assert "theta should be between [0, pi]" in str(exc_info.value) 653 654 with assert_raises(ValueError) as exc_info: 655 invalid_theta = linspace(0.1, 1.1, num=ntheta) * pi 656 invalid_lats, lons = meshgrid(invalid_theta, phi) 657 LSQSphereBivariateSpline(invalid_lats.ravel(), lons.ravel(), 658 data.T.ravel(), knotst, knotsp) 659 assert "theta should be between [0, pi]" in str(exc_info.value) 660 661 with assert_raises(ValueError) as exc_info: 662 invalid_phi = linspace(-0.1, 1.0, num=ntheta) * 2.0 * pi 663 lats, invalid_lons = meshgrid(theta, invalid_phi) 664 LSQSphereBivariateSpline(lats.ravel(), invalid_lons.ravel(), 665 data.T.ravel(), knotst, knotsp) 666 assert "phi should be between [0, 2pi]" in str(exc_info.value) 667 668 with assert_raises(ValueError) as exc_info: 669 invalid_phi = linspace(0.0, 1.1, num=ntheta) * 2.0 * pi 670 lats, invalid_lons = meshgrid(theta, invalid_phi) 671 LSQSphereBivariateSpline(lats.ravel(), invalid_lons.ravel(), 672 data.T.ravel(), knotst, knotsp) 673 assert "phi should be between [0, 2pi]" in str(exc_info.value) 674 675 lats, lons = meshgrid(theta, phi) 676 677 with assert_raises(ValueError) as exc_info: 678 invalid_knotst = np.copy(knotst) 679 invalid_knotst[0] = -0.1 680 LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), 681 data.T.ravel(), invalid_knotst, knotsp) 682 assert "tt should be between (0, pi)" in str(exc_info.value) 683 684 with assert_raises(ValueError) as exc_info: 685 invalid_knotst = np.copy(knotst) 686 invalid_knotst[0] = pi 687 LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), 688 data.T.ravel(), invalid_knotst, knotsp) 689 assert "tt should be between (0, pi)" in str(exc_info.value) 690 691 with assert_raises(ValueError) as exc_info: 692 invalid_knotsp = np.copy(knotsp) 693 invalid_knotsp[0] = -0.1 694 LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), 695 data.T.ravel(), knotst, invalid_knotsp) 696 assert "tp should be between (0, 2pi)" in str(exc_info.value) 697 698 with assert_raises(ValueError) as exc_info: 699 invalid_knotsp = np.copy(knotsp) 700 invalid_knotsp[0] = 2 * pi 701 LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), 702 data.T.ravel(), knotst, invalid_knotsp) 703 assert "tp should be between (0, 2pi)" in str(exc_info.value) 704 705 with assert_raises(ValueError) as exc_info: 706 invalid_w = array([-1.0, 1.0, 1.5, 0.5, 1.0, 1.5, 0.5, 1.0, 1.0]) 707 LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), data.T.ravel(), 708 knotst, knotsp, w=invalid_w) 709 assert "w should be positive" in str(exc_info.value) 710 711 with assert_raises(ValueError) as exc_info: 712 LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), data.T.ravel(), 713 knotst, knotsp, eps=0.0) 714 assert "eps should be between (0, 1)" in str(exc_info.value) 715 716 with assert_raises(ValueError) as exc_info: 717 LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), data.T.ravel(), 718 knotst, knotsp, eps=1.0) 719 assert "eps should be between (0, 1)" in str(exc_info.value) 720 721 def test_array_like_input(self): 722 ntheta, nphi = 70, 90 723 theta = linspace(0.5 / (ntheta - 1), 1 - 0.5 / (ntheta - 1), 724 ntheta) * pi 725 phi = linspace(0.5 / (nphi - 1), 1 - 0.5 / (nphi - 1), 726 nphi) * 2. * pi 727 lats, lons = meshgrid(theta, phi) 728 data = ones((theta.shape[0], phi.shape[0])) 729 # define knots and extract data values at the knots 730 knotst = theta[::5] 731 knotsp = phi[::5] 732 w = ones((lats.ravel().shape[0])) 733 734 # np.array input 735 spl1 = LSQSphereBivariateSpline(lats.ravel(), lons.ravel(), 736 data.T.ravel(), knotst, knotsp, w=w) 737 # list input 738 spl2 = LSQSphereBivariateSpline(lats.ravel().tolist(), 739 lons.ravel().tolist(), 740 data.T.ravel().tolist(), 741 knotst.tolist(), 742 knotsp.tolist(), w=w.tolist()) 743 assert_array_almost_equal(spl1(1.0, 1.0), spl2(1.0, 1.0)) 744 745 746class TestSmoothSphereBivariateSpline: 747 def setup_method(self): 748 theta = array([.25*pi, .25*pi, .25*pi, .5*pi, .5*pi, .5*pi, .75*pi, 749 .75*pi, .75*pi]) 750 phi = array([.5 * pi, pi, 1.5 * pi, .5 * pi, pi, 1.5 * pi, .5 * pi, pi, 751 1.5 * pi]) 752 r = array([3, 3, 3, 3, 3, 3, 3, 3, 3]) 753 self.lut = SmoothSphereBivariateSpline(theta, phi, r, s=1E10) 754 755 def test_linear_constant(self): 756 assert_almost_equal(self.lut.get_residual(), 0.) 757 assert_array_almost_equal(self.lut([1, 1.5, 2],[1, 1.5]), 758 [[3, 3], [3, 3], [3, 3]]) 759 760 def test_empty_input(self): 761 assert_array_almost_equal(self.lut([], []), np.zeros((0,0))) 762 assert_array_almost_equal(self.lut([], [], grid=False), np.zeros((0,))) 763 764 def test_invalid_input(self): 765 theta = array([.25 * pi, .25 * pi, .25 * pi, .5 * pi, .5 * pi, .5 * pi, 766 .75 * pi, .75 * pi, .75 * pi]) 767 phi = array([.5 * pi, pi, 1.5 * pi, .5 * pi, pi, 1.5 * pi, .5 * pi, pi, 768 1.5 * pi]) 769 r = array([3, 3, 3, 3, 3, 3, 3, 3, 3]) 770 771 with assert_raises(ValueError) as exc_info: 772 invalid_theta = array([-0.1 * pi, .25 * pi, .25 * pi, .5 * pi, 773 .5 * pi, .5 * pi, .75 * pi, .75 * pi, 774 .75 * pi]) 775 SmoothSphereBivariateSpline(invalid_theta, phi, r, s=1E10) 776 assert "theta should be between [0, pi]" in str(exc_info.value) 777 778 with assert_raises(ValueError) as exc_info: 779 invalid_theta = array([.25 * pi, .25 * pi, .25 * pi, .5 * pi, 780 .5 * pi, .5 * pi, .75 * pi, .75 * pi, 781 1.1 * pi]) 782 SmoothSphereBivariateSpline(invalid_theta, phi, r, s=1E10) 783 assert "theta should be between [0, pi]" in str(exc_info.value) 784 785 with assert_raises(ValueError) as exc_info: 786 invalid_phi = array([-.1 * pi, pi, 1.5 * pi, .5 * pi, pi, 1.5 * pi, 787 .5 * pi, pi, 1.5 * pi]) 788 SmoothSphereBivariateSpline(theta, invalid_phi, r, s=1E10) 789 assert "phi should be between [0, 2pi]" in str(exc_info.value) 790 791 with assert_raises(ValueError) as exc_info: 792 invalid_phi = array([1.0 * pi, pi, 1.5 * pi, .5 * pi, pi, 1.5 * pi, 793 .5 * pi, pi, 2.1 * pi]) 794 SmoothSphereBivariateSpline(theta, invalid_phi, r, s=1E10) 795 assert "phi should be between [0, 2pi]" in str(exc_info.value) 796 797 with assert_raises(ValueError) as exc_info: 798 invalid_w = array([-1.0, 1.0, 1.5, 0.5, 1.0, 1.5, 0.5, 1.0, 1.0]) 799 SmoothSphereBivariateSpline(theta, phi, r, w=invalid_w, s=1E10) 800 assert "w should be positive" in str(exc_info.value) 801 802 with assert_raises(ValueError) as exc_info: 803 SmoothSphereBivariateSpline(theta, phi, r, s=-1.0) 804 assert "s should be positive" in str(exc_info.value) 805 806 with assert_raises(ValueError) as exc_info: 807 SmoothSphereBivariateSpline(theta, phi, r, eps=-1.0) 808 assert "eps should be between (0, 1)" in str(exc_info.value) 809 810 with assert_raises(ValueError) as exc_info: 811 SmoothSphereBivariateSpline(theta, phi, r, eps=1.0) 812 assert "eps should be between (0, 1)" in str(exc_info.value) 813 814 def test_array_like_input(self): 815 theta = np.array([.25 * pi, .25 * pi, .25 * pi, .5 * pi, .5 * pi, 816 .5 * pi, .75 * pi, .75 * pi, .75 * pi]) 817 phi = np.array([.5 * pi, pi, 1.5 * pi, .5 * pi, pi, 1.5 * pi, .5 * pi, 818 pi, 1.5 * pi]) 819 r = np.array([3, 3, 3, 3, 3, 3, 3, 3, 3]) 820 w = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) 821 822 # np.array input 823 spl1 = SmoothSphereBivariateSpline(theta, phi, r, w=w, s=1E10) 824 825 # list input 826 spl2 = SmoothSphereBivariateSpline(theta.tolist(), phi.tolist(), 827 r.tolist(), w=w.tolist(), s=1E10) 828 assert_array_almost_equal(spl1(1.0, 1.0), spl2(1.0, 1.0)) 829 830 831class TestRectBivariateSpline: 832 def test_defaults(self): 833 x = array([1,2,3,4,5]) 834 y = array([1,2,3,4,5]) 835 z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]]) 836 lut = RectBivariateSpline(x,y,z) 837 assert_array_almost_equal(lut(x,y),z) 838 839 def test_evaluate(self): 840 x = array([1,2,3,4,5]) 841 y = array([1,2,3,4,5]) 842 z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]]) 843 lut = RectBivariateSpline(x,y,z) 844 845 xi = [1, 2.3, 5.3, 0.5, 3.3, 1.2, 3] 846 yi = [1, 3.3, 1.2, 4.0, 5.0, 1.0, 3] 847 zi = lut.ev(xi, yi) 848 zi2 = array([lut(xp, yp)[0,0] for xp, yp in zip(xi, yi)]) 849 850 assert_almost_equal(zi, zi2) 851 852 def test_derivatives_grid(self): 853 x = array([1,2,3,4,5]) 854 y = array([1,2,3,4,5]) 855 z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]]) 856 dx = array([[0,0,-20,0,0],[0,0,13,0,0],[0,0,4,0,0], 857 [0,0,-11,0,0],[0,0,4,0,0]])/6. 858 dy = array([[4,-1,0,1,-4],[4,-1,0,1,-4],[0,1.5,0,-1.5,0], 859 [2,.25,0,-.25,-2],[4,-1,0,1,-4]]) 860 dxdy = array([[40,-25,0,25,-40],[-26,16.25,0,-16.25,26], 861 [-8,5,0,-5,8],[22,-13.75,0,13.75,-22],[-8,5,0,-5,8]])/6. 862 lut = RectBivariateSpline(x,y,z) 863 assert_array_almost_equal(lut(x,y,dx=1),dx) 864 assert_array_almost_equal(lut(x,y,dy=1),dy) 865 assert_array_almost_equal(lut(x,y,dx=1,dy=1),dxdy) 866 867 def test_derivatives(self): 868 x = array([1,2,3,4,5]) 869 y = array([1,2,3,4,5]) 870 z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]]) 871 dx = array([0,0,2./3,0,0]) 872 dy = array([4,-1,0,-.25,-4]) 873 dxdy = array([160,65,0,55,32])/24. 874 lut = RectBivariateSpline(x,y,z) 875 assert_array_almost_equal(lut(x,y,dx=1,grid=False),dx) 876 assert_array_almost_equal(lut(x,y,dy=1,grid=False),dy) 877 assert_array_almost_equal(lut(x,y,dx=1,dy=1,grid=False),dxdy) 878 879 def test_broadcast(self): 880 x = array([1,2,3,4,5]) 881 y = array([1,2,3,4,5]) 882 z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]]) 883 lut = RectBivariateSpline(x,y,z) 884 assert_allclose(lut(x, y), lut(x[:,None], y[None,:], grid=False)) 885 886 def test_invalid_input(self): 887 888 with assert_raises(ValueError) as info: 889 x = array([6, 2, 3, 4, 5]) 890 y = array([1, 2, 3, 4, 5]) 891 z = array([[1, 2, 1, 2, 1], [1, 2, 1, 2, 1], [1, 2, 3, 2, 1], 892 [1, 2, 2, 2, 1], [1, 2, 1, 2, 1]]) 893 RectBivariateSpline(x, y, z) 894 assert "x must be strictly increasing" in str(info.value) 895 896 with assert_raises(ValueError) as info: 897 x = array([1, 2, 3, 4, 5]) 898 y = array([2, 2, 3, 4, 5]) 899 z = array([[1, 2, 1, 2, 1], [1, 2, 1, 2, 1], [1, 2, 3, 2, 1], 900 [1, 2, 2, 2, 1], [1, 2, 1, 2, 1]]) 901 RectBivariateSpline(x, y, z) 902 assert "y must be strictly increasing" in str(info.value) 903 904 with assert_raises(ValueError) as info: 905 x = array([1, 2, 3, 4, 5]) 906 y = array([1, 2, 3, 4, 5]) 907 z = array([[1, 2, 1, 2, 1], [1, 2, 1, 2, 1], [1, 2, 3, 2, 1], 908 [1, 2, 2, 2, 1]]) 909 RectBivariateSpline(x, y, z) 910 assert "x dimension of z must have same number of elements as x"\ 911 in str(info.value) 912 913 with assert_raises(ValueError) as info: 914 x = array([1, 2, 3, 4, 5]) 915 y = array([1, 2, 3, 4, 5]) 916 z = array([[1, 2, 1, 2], [1, 2, 1, 2], [1, 2, 3, 2], 917 [1, 2, 2, 2], [1, 2, 1, 2]]) 918 RectBivariateSpline(x, y, z) 919 assert "y dimension of z must have same number of elements as y"\ 920 in str(info.value) 921 922 with assert_raises(ValueError) as info: 923 x = array([1, 2, 3, 4, 5]) 924 y = array([1, 2, 3, 4, 5]) 925 z = array([[1, 2, 1, 2, 1], [1, 2, 1, 2, 1], [1, 2, 3, 2, 1], 926 [1, 2, 2, 2, 1], [1, 2, 1, 2, 1]]) 927 bbox = (-100, 100, -100) 928 RectBivariateSpline(x, y, z, bbox=bbox) 929 assert "bbox shape should be (4,)" in str(info.value) 930 931 with assert_raises(ValueError) as info: 932 RectBivariateSpline(x, y, z, s=-1.0) 933 assert "s should be s >= 0.0" in str(info.value) 934 935 def test_array_like_input(self): 936 x = array([1, 2, 3, 4, 5]) 937 y = array([1, 2, 3, 4, 5]) 938 z = array([[1, 2, 1, 2, 1], [1, 2, 1, 2, 1], [1, 2, 3, 2, 1], 939 [1, 2, 2, 2, 1], [1, 2, 1, 2, 1]]) 940 bbox = array([1, 5, 1, 5]) 941 942 spl1 = RectBivariateSpline(x, y, z, bbox=bbox) 943 spl2 = RectBivariateSpline(x.tolist(), y.tolist(), z.tolist(), 944 bbox=bbox.tolist()) 945 assert_array_almost_equal(spl1(1.0, 1.0), spl2(1.0, 1.0)) 946 947 def test_not_increasing_input(self): 948 # gh-8565 949 NSamp = 20 950 Theta = np.random.uniform(0, np.pi, NSamp) 951 Phi = np.random.uniform(0, 2 * np.pi, NSamp) 952 Data = np.ones(NSamp) 953 954 Interpolator = SmoothSphereBivariateSpline(Theta, Phi, Data, s=3.5) 955 956 NLon = 6 957 NLat = 3 958 GridPosLats = np.arange(NLat) / NLat * np.pi 959 GridPosLons = np.arange(NLon) / NLon * 2 * np.pi 960 961 # No error 962 Interpolator(GridPosLats, GridPosLons) 963 964 nonGridPosLats = GridPosLats.copy() 965 nonGridPosLats[2] = 0.001 966 with assert_raises(ValueError) as exc_info: 967 Interpolator(nonGridPosLats, GridPosLons) 968 assert "x must be strictly increasing" in str(exc_info.value) 969 970 nonGridPosLons = GridPosLons.copy() 971 nonGridPosLons[2] = 0.001 972 with assert_raises(ValueError) as exc_info: 973 Interpolator(GridPosLats, nonGridPosLons) 974 assert "y must be strictly increasing" in str(exc_info.value) 975 976 977class TestRectSphereBivariateSpline: 978 def test_defaults(self): 979 y = linspace(0.01, 2*pi-0.01, 7) 980 x = linspace(0.01, pi-0.01, 7) 981 z = array([[1,2,1,2,1,2,1],[1,2,1,2,1,2,1],[1,2,3,2,1,2,1], 982 [1,2,2,2,1,2,1],[1,2,1,2,1,2,1],[1,2,2,2,1,2,1], 983 [1,2,1,2,1,2,1]]) 984 lut = RectSphereBivariateSpline(x,y,z) 985 assert_array_almost_equal(lut(x,y),z) 986 987 def test_evaluate(self): 988 y = linspace(0.01, 2*pi-0.01, 7) 989 x = linspace(0.01, pi-0.01, 7) 990 z = array([[1,2,1,2,1,2,1],[1,2,1,2,1,2,1],[1,2,3,2,1,2,1], 991 [1,2,2,2,1,2,1],[1,2,1,2,1,2,1],[1,2,2,2,1,2,1], 992 [1,2,1,2,1,2,1]]) 993 lut = RectSphereBivariateSpline(x,y,z) 994 yi = [0.2, 1, 2.3, 2.35, 3.0, 3.99, 5.25] 995 xi = [1.5, 0.4, 1.1, 0.45, 0.2345, 1., 0.0001] 996 zi = lut.ev(xi, yi) 997 zi2 = array([lut(xp, yp)[0,0] for xp, yp in zip(xi, yi)]) 998 assert_almost_equal(zi, zi2) 999 1000 def test_derivatives_grid(self): 1001 y = linspace(0.01, 2*pi-0.01, 7) 1002 x = linspace(0.01, pi-0.01, 7) 1003 z = array([[1,2,1,2,1,2,1],[1,2,1,2,1,2,1],[1,2,3,2,1,2,1], 1004 [1,2,2,2,1,2,1],[1,2,1,2,1,2,1],[1,2,2,2,1,2,1], 1005 [1,2,1,2,1,2,1]]) 1006 1007 lut = RectSphereBivariateSpline(x,y,z) 1008 1009 y = linspace(0.02, 2*pi-0.02, 7) 1010 x = linspace(0.02, pi-0.02, 7) 1011 1012 assert_allclose(lut(x, y, dtheta=1), _numdiff_2d(lut, x, y, dx=1), 1013 rtol=1e-4, atol=1e-4) 1014 assert_allclose(lut(x, y, dphi=1), _numdiff_2d(lut, x, y, dy=1), 1015 rtol=1e-4, atol=1e-4) 1016 assert_allclose(lut(x, y, dtheta=1, dphi=1), _numdiff_2d(lut, x, y, dx=1, dy=1, eps=1e-6), 1017 rtol=1e-3, atol=1e-3) 1018 1019 def test_derivatives(self): 1020 y = linspace(0.01, 2*pi-0.01, 7) 1021 x = linspace(0.01, pi-0.01, 7) 1022 z = array([[1,2,1,2,1,2,1],[1,2,1,2,1,2,1],[1,2,3,2,1,2,1], 1023 [1,2,2,2,1,2,1],[1,2,1,2,1,2,1],[1,2,2,2,1,2,1], 1024 [1,2,1,2,1,2,1]]) 1025 1026 lut = RectSphereBivariateSpline(x,y,z) 1027 1028 y = linspace(0.02, 2*pi-0.02, 7) 1029 x = linspace(0.02, pi-0.02, 7) 1030 1031 assert_equal(lut(x, y, dtheta=1, grid=False).shape, x.shape) 1032 assert_allclose(lut(x, y, dtheta=1, grid=False), 1033 _numdiff_2d(lambda x,y: lut(x,y,grid=False), x, y, dx=1), 1034 rtol=1e-4, atol=1e-4) 1035 assert_allclose(lut(x, y, dphi=1, grid=False), 1036 _numdiff_2d(lambda x,y: lut(x,y,grid=False), x, y, dy=1), 1037 rtol=1e-4, atol=1e-4) 1038 assert_allclose(lut(x, y, dtheta=1, dphi=1, grid=False), 1039 _numdiff_2d(lambda x,y: lut(x,y,grid=False), x, y, dx=1, dy=1, eps=1e-6), 1040 rtol=1e-3, atol=1e-3) 1041 1042 def test_invalid_input(self): 1043 data = np.dot(np.atleast_2d(90. - np.linspace(-80., 80., 18)).T, 1044 np.atleast_2d(180. - np.abs(np.linspace(0., 350., 9)))).T 1045 1046 with assert_raises(ValueError) as exc_info: 1047 lats = np.linspace(0, 170, 9) * np.pi / 180. 1048 lons = np.linspace(0, 350, 18) * np.pi / 180. 1049 RectSphereBivariateSpline(lats, lons, data) 1050 assert "u should be between (0, pi)" in str(exc_info.value) 1051 1052 with assert_raises(ValueError) as exc_info: 1053 lats = np.linspace(10, 180, 9) * np.pi / 180. 1054 lons = np.linspace(0, 350, 18) * np.pi / 180. 1055 RectSphereBivariateSpline(lats, lons, data) 1056 assert "u should be between (0, pi)" in str(exc_info.value) 1057 1058 with assert_raises(ValueError) as exc_info: 1059 lats = np.linspace(10, 170, 9) * np.pi / 180. 1060 lons = np.linspace(-181, 10, 18) * np.pi / 180. 1061 RectSphereBivariateSpline(lats, lons, data) 1062 assert "v[0] should be between [-pi, pi)" in str(exc_info.value) 1063 1064 with assert_raises(ValueError) as exc_info: 1065 lats = np.linspace(10, 170, 9) * np.pi / 180. 1066 lons = np.linspace(-10, 360, 18) * np.pi / 180. 1067 RectSphereBivariateSpline(lats, lons, data) 1068 assert "v[-1] should be v[0] + 2pi or less" in str(exc_info.value) 1069 1070 with assert_raises(ValueError) as exc_info: 1071 lats = np.linspace(10, 170, 9) * np.pi / 180. 1072 lons = np.linspace(10, 350, 18) * np.pi / 180. 1073 RectSphereBivariateSpline(lats, lons, data, s=-1) 1074 assert "s should be positive" in str(exc_info.value) 1075 1076 def test_array_like_input(self): 1077 y = linspace(0.01, 2 * pi - 0.01, 7) 1078 x = linspace(0.01, pi - 0.01, 7) 1079 z = array([[1, 2, 1, 2, 1, 2, 1], [1, 2, 1, 2, 1, 2, 1], 1080 [1, 2, 3, 2, 1, 2, 1], 1081 [1, 2, 2, 2, 1, 2, 1], [1, 2, 1, 2, 1, 2, 1], 1082 [1, 2, 2, 2, 1, 2, 1], 1083 [1, 2, 1, 2, 1, 2, 1]]) 1084 # np.array input 1085 spl1 = RectSphereBivariateSpline(x, y, z) 1086 # list input 1087 spl2 = RectSphereBivariateSpline(x.tolist(), y.tolist(), z.tolist()) 1088 assert_array_almost_equal(spl1(x, y), spl2(x, y)) 1089 1090 def test_negative_evaluation(self): 1091 lats = np.array([25, 30, 35, 40, 45]) 1092 lons = np.array([-90, -85, -80, -75, 70]) 1093 mesh = np.meshgrid(lats, lons) 1094 data = mesh[0] + mesh[1] # lon + lat value 1095 lat_r = np.radians(lats) 1096 lon_r = np.radians(lons) 1097 interpolator = RectSphereBivariateSpline(lat_r, lon_r, data) 1098 query_lat = np.radians(np.array([35, 37.5])) 1099 query_lon = np.radians(np.array([-80, -77.5])) 1100 data_interp = interpolator(query_lat, query_lon) 1101 ans = np.array([[-45.0, -42.480862], 1102 [-49.0625, -46.54315]]) 1103 assert_array_almost_equal(data_interp, ans) 1104 1105 1106def _numdiff_2d(func, x, y, dx=0, dy=0, eps=1e-8): 1107 if dx == 0 and dy == 0: 1108 return func(x, y) 1109 elif dx == 1 and dy == 0: 1110 return (func(x + eps, y) - func(x - eps, y)) / (2*eps) 1111 elif dx == 0 and dy == 1: 1112 return (func(x, y + eps) - func(x, y - eps)) / (2*eps) 1113 elif dx == 1 and dy == 1: 1114 return (func(x + eps, y + eps) - func(x - eps, y + eps) 1115 - func(x + eps, y - eps) + func(x - eps, y - eps)) / (2*eps)**2 1116 else: 1117 raise ValueError("invalid derivative order") 1118