1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2 3import pytest 4import numpy as np 5import numpy.ma as ma 6from contextlib import nullcontext 7 8from astropy.convolution.convolve import convolve, convolve_fft 9from astropy.convolution.kernels import Gaussian2DKernel 10from astropy.utils.exceptions import AstropyUserWarning 11from astropy import units as u 12from astropy.utils.compat.optional_deps import HAS_SCIPY, HAS_PANDAS # noqa 13 14from numpy.testing import (assert_array_almost_equal_nulp, 15 assert_array_almost_equal, 16 assert_allclose) 17 18import itertools 19 20VALID_DTYPES = ('>f4', '<f4', '>f8', '<f8') 21VALID_DTYPE_MATRIX = list(itertools.product(VALID_DTYPES, VALID_DTYPES)) 22 23BOUNDARY_OPTIONS = [None, 'fill', 'wrap', 'extend'] 24NANHANDLING_OPTIONS = ['interpolate', 'fill'] 25NORMALIZE_OPTIONS = [True, False] 26PRESERVE_NAN_OPTIONS = [True, False] 27 28BOUNDARIES_AND_CONVOLUTIONS = (list(zip(itertools.cycle((convolve,)), 29 BOUNDARY_OPTIONS)) + [(convolve_fft, 30 'wrap'), 31 (convolve_fft, 32 'fill')]) 33 34 35class TestConvolve1D: 36 37 def test_list(self): 38 """ 39 Test that convolve works correctly when inputs are lists 40 """ 41 42 x = [1, 4, 5, 6, 5, 7, 8] 43 y = [0.2, 0.6, 0.2] 44 z = convolve(x, y, boundary=None) 45 assert_array_almost_equal_nulp(z, 46 np.array([0., 3.6, 5., 5.6, 5.6, 6.8, 0.]), 10) 47 48 def test_tuple(self): 49 """ 50 Test that convolve works correctly when inputs are tuples 51 """ 52 53 x = (1, 4, 5, 6, 5, 7, 8) 54 y = (0.2, 0.6, 0.2) 55 z = convolve(x, y, boundary=None) 56 assert_array_almost_equal_nulp(z, 57 np.array([0., 3.6, 5., 5.6, 5.6, 6.8, 0.]), 10) 58 59 @pytest.mark.parametrize(('boundary', 'nan_treatment', 60 'normalize_kernel', 'preserve_nan', 'dtype'), 61 itertools.product(BOUNDARY_OPTIONS, 62 NANHANDLING_OPTIONS, 63 NORMALIZE_OPTIONS, 64 PRESERVE_NAN_OPTIONS, 65 VALID_DTYPES)) 66 def test_quantity(self, boundary, nan_treatment, 67 normalize_kernel, preserve_nan, dtype): 68 """ 69 Test that convolve works correctly when input array is a Quantity 70 """ 71 72 x = np.array([1, 4, 5, 6, 5, 7, 8], dtype=dtype) * u.ph 73 y = np.array([0.2, 0.6, 0.2], dtype=dtype) 74 z = convolve(x, y, boundary=boundary, nan_treatment=nan_treatment, 75 normalize_kernel=normalize_kernel, preserve_nan=preserve_nan) 76 77 assert x.unit == z.unit 78 79 @pytest.mark.parametrize(('boundary', 'nan_treatment', 80 'normalize_kernel', 'preserve_nan', 'dtype'), 81 itertools.product(BOUNDARY_OPTIONS, 82 NANHANDLING_OPTIONS, 83 NORMALIZE_OPTIONS, 84 PRESERVE_NAN_OPTIONS, 85 VALID_DTYPES)) 86 def test_input_unmodified(self, boundary, nan_treatment, 87 normalize_kernel, preserve_nan, dtype): 88 """ 89 Test that convolve works correctly when inputs are lists 90 """ 91 92 array = [1., 4., 5., 6., 5., 7., 8.] 93 kernel = [0.2, 0.6, 0.2] 94 x = np.array(array, dtype=dtype) 95 y = np.array(kernel, dtype=dtype) 96 97 # Make pseudoimmutable 98 x.flags.writeable = False 99 y.flags.writeable = False 100 101 z = convolve(x, y, boundary=boundary, nan_treatment=nan_treatment, 102 normalize_kernel=normalize_kernel, preserve_nan=preserve_nan) 103 104 assert np.all(np.array(array, dtype=dtype) == x) 105 assert np.all(np.array(kernel, dtype=dtype) == y) 106 107 @pytest.mark.parametrize(('boundary', 'nan_treatment', 108 'normalize_kernel', 'preserve_nan', 'dtype'), 109 itertools.product(BOUNDARY_OPTIONS, 110 NANHANDLING_OPTIONS, 111 NORMALIZE_OPTIONS, 112 PRESERVE_NAN_OPTIONS, 113 VALID_DTYPES)) 114 def test_input_unmodified_with_nan(self, boundary, nan_treatment, 115 normalize_kernel, preserve_nan, dtype): 116 """ 117 Test that convolve doesn't modify the input data 118 """ 119 120 array = [1., 4., 5., np.nan, 5., 7., 8.] 121 kernel = [0.2, 0.6, 0.2] 122 x = np.array(array, dtype=dtype) 123 y = np.array(kernel, dtype=dtype) 124 125 # Make pseudoimmutable 126 x.flags.writeable = False 127 y.flags.writeable = False 128 129 # make copies for post call comparison 130 x_copy = x.copy() 131 y_copy = y.copy() 132 133 z = convolve(x, y, boundary=boundary, nan_treatment=nan_treatment, 134 normalize_kernel=normalize_kernel, preserve_nan=preserve_nan) 135 136 # ( NaN == NaN ) = False 137 # Only compare non NaN values for canonical equivalence 138 # and then check NaN explicitly with np.isnan() 139 array_is_nan = np.isnan(array) 140 kernel_is_nan = np.isnan(kernel) 141 array_not_nan = ~array_is_nan 142 kernel_not_nan = ~kernel_is_nan 143 assert np.all(x_copy[array_not_nan] == x[array_not_nan]) 144 assert np.all(y_copy[kernel_not_nan] == y[kernel_not_nan]) 145 assert np.all(np.isnan(x[array_is_nan])) 146 assert np.all(np.isnan(y[kernel_is_nan])) 147 148 @pytest.mark.parametrize(('dtype_array', 'dtype_kernel'), VALID_DTYPE_MATRIX) 149 def test_dtype(self, dtype_array, dtype_kernel): 150 ''' 151 Test that 32- and 64-bit floats are correctly handled 152 ''' 153 154 x = np.array([1., 2., 3.], dtype=dtype_array) 155 156 y = np.array([0., 1., 0.], dtype=dtype_kernel) 157 158 z = convolve(x, y) 159 160 assert x.dtype == z.dtype 161 162 @pytest.mark.parametrize(('convfunc', 'boundary',), BOUNDARIES_AND_CONVOLUTIONS) 163 def test_unity_1_none(self, boundary, convfunc): 164 ''' 165 Test that a unit kernel with a single element returns the same array 166 ''' 167 168 x = np.array([1., 2., 3.], dtype='>f8') 169 170 y = np.array([1.], dtype='>f8') 171 172 z = convfunc(x, y, boundary=boundary) 173 174 np.testing.assert_allclose(z, x) 175 176 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 177 def test_unity_3(self, boundary): 178 ''' 179 Test that a unit kernel with three elements returns the same array 180 (except when boundary is None). 181 ''' 182 183 x = np.array([1., 2., 3.], dtype='>f8') 184 185 y = np.array([0., 1., 0.], dtype='>f8') 186 187 z = convolve(x, y, boundary=boundary) 188 189 if boundary is None: 190 assert np.all(z == np.array([0., 2., 0.], dtype='>f8')) 191 else: 192 assert np.all(z == x) 193 194 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 195 def test_uniform_3(self, boundary): 196 ''' 197 Test that the different modes are producing the correct results using 198 a uniform kernel with three elements 199 ''' 200 201 x = np.array([1., 0., 3.], dtype='>f8') 202 203 y = np.array([1., 1., 1.], dtype='>f8') 204 205 z = convolve(x, y, boundary=boundary, normalize_kernel=False) 206 207 if boundary is None: 208 assert np.all(z == np.array([0., 4., 0.], dtype='>f8')) 209 elif boundary == 'fill': 210 assert np.all(z == np.array([1., 4., 3.], dtype='>f8')) 211 elif boundary == 'wrap': 212 assert np.all(z == np.array([4., 4., 4.], dtype='>f8')) 213 else: 214 assert np.all(z == np.array([2., 4., 6.], dtype='>f8')) 215 216 @pytest.mark.parametrize(('boundary', 'nan_treatment', 217 'normalize_kernel', 'preserve_nan'), 218 itertools.product(BOUNDARY_OPTIONS, 219 NANHANDLING_OPTIONS, 220 NORMALIZE_OPTIONS, 221 PRESERVE_NAN_OPTIONS)) 222 def test_unity_3_withnan(self, boundary, nan_treatment, 223 normalize_kernel, preserve_nan): 224 ''' 225 Test that a unit kernel with three elements returns the same array 226 (except when boundary is None). This version includes a NaN value in 227 the original array. 228 ''' 229 230 x = np.array([1., np.nan, 3.], dtype='>f8') 231 232 y = np.array([0., 1., 0.], dtype='>f8') 233 234 # Since the kernel is actually only one pixel wide (because of the 235 # zeros) the NaN value doesn't get interpolated over so a warning is 236 # expected. 237 if nan_treatment == 'interpolate' and not preserve_nan: 238 ctx = pytest.warns(AstropyUserWarning, 239 match="nan_treatment='interpolate', however, " 240 "NaN values detected") 241 else: 242 ctx = nullcontext() 243 244 with ctx: 245 z = convolve(x, y, boundary=boundary, nan_treatment=nan_treatment, 246 normalize_kernel=normalize_kernel, 247 preserve_nan=preserve_nan) 248 249 if preserve_nan: 250 assert np.isnan(z[1]) 251 252 x = np.nan_to_num(z) 253 z = np.nan_to_num(z) 254 255 if boundary is None: 256 assert np.all(z == np.array([0., 0., 0.], dtype='>f8')) 257 else: 258 assert np.all(z == x) 259 260 @pytest.mark.parametrize(('boundary', 'nan_treatment', 261 'normalize_kernel', 'preserve_nan'), 262 itertools.product(BOUNDARY_OPTIONS, 263 NANHANDLING_OPTIONS, 264 NORMALIZE_OPTIONS, 265 PRESERVE_NAN_OPTIONS)) 266 def test_uniform_3_withnan(self, boundary, nan_treatment, normalize_kernel, 267 preserve_nan): 268 ''' 269 Test that the different modes are producing the correct results using 270 a uniform kernel with three elements. This version includes a NaN 271 value in the original array. 272 ''' 273 274 x = np.array([1., np.nan, 3.], dtype='>f8') 275 276 y = np.array([1., 1., 1.], dtype='>f8') 277 278 z = convolve(x, y, boundary=boundary, nan_treatment=nan_treatment, 279 normalize_kernel=normalize_kernel, 280 preserve_nan=preserve_nan) 281 282 if preserve_nan: 283 assert np.isnan(z[1]) 284 285 z = np.nan_to_num(z) 286 287 # boundary, nan_treatment, normalize_kernel 288 rslt = { 289 (None, 'interpolate', True): [0, 2, 0], 290 (None, 'interpolate', False): [0, 6, 0], 291 (None, 'fill', True): [0, 4/3., 0], 292 (None, 'fill', False): [0, 4, 0], 293 ('fill', 'interpolate', True): [1/2., 2, 3/2.], 294 ('fill', 'interpolate', False): [3/2., 6, 9/2.], 295 ('fill', 'fill', True): [1/3., 4/3., 3/3.], 296 ('fill', 'fill', False): [1, 4, 3], 297 ('wrap', 'interpolate', True): [2, 2, 2], 298 ('wrap', 'interpolate', False): [6, 6, 6], 299 ('wrap', 'fill', True): [4/3., 4/3., 4/3.], 300 ('wrap', 'fill', False): [4, 4, 4], 301 ('extend', 'interpolate', True): [1, 2, 3], 302 ('extend', 'interpolate', False): [3, 6, 9], 303 ('extend', 'fill', True): [2/3., 4/3., 6/3.], 304 ('extend', 'fill', False): [2, 4, 6], 305 }[boundary, nan_treatment, normalize_kernel] 306 if preserve_nan: 307 rslt[1] = 0 308 309 assert_array_almost_equal_nulp(z, np.array(rslt, dtype='>f8'), 10) 310 311 @pytest.mark.parametrize(('boundary', 'normalize_kernel'), 312 itertools.product(BOUNDARY_OPTIONS, 313 NORMALIZE_OPTIONS)) 314 def test_zero_sum_kernel(self, boundary, normalize_kernel): 315 """ 316 Test that convolve works correctly with zero sum kernels. 317 """ 318 319 if normalize_kernel: 320 pytest.xfail("You can't normalize by a zero sum kernel") 321 322 x = [1, 2, 3, 4, 5, 6, 7, 8, 9] 323 y = [-1, -1, -1, -1, 8, -1, -1, -1, -1] 324 assert(np.isclose(sum(y), 0, atol=1e-8)) 325 326 z = convolve(x, y, boundary=boundary, normalize_kernel=normalize_kernel) 327 328 # boundary, normalize_kernel == False 329 rslt = { 330 (None): [0., 0., 0., 0., 0., 0., 0., 0., 0.], 331 ('fill'): [-6., -3., -1., 0., 0., 10., 21., 33., 46.], 332 ('wrap'): [-36., -27., -18., -9., 0., 9., 18., 27., 36.], 333 ('extend'): [-10., -6., -3., -1., 0., 1., 3., 6., 10.] 334 }[boundary] 335 336 assert_array_almost_equal_nulp(z, np.array(rslt, dtype='>f8'), 10) 337 338 @pytest.mark.parametrize(('boundary', 'normalize_kernel'), 339 itertools.product(BOUNDARY_OPTIONS, 340 NORMALIZE_OPTIONS)) 341 def test_int_masked_kernel(self, boundary, normalize_kernel): 342 """ 343 Test that convolve works correctly with integer masked kernels. 344 """ 345 346 if normalize_kernel: 347 pytest.xfail("You can't normalize by a zero sum kernel") 348 349 x = [1, 2, 3, 4, 5, 6, 7, 8, 9] 350 y = ma.array([-1, -1, -1, -1, 8, -1, -1, -1, -1], mask=[1, 0, 0, 0, 0, 0, 0, 0, 0], fill_value=0.) 351 352 z = convolve(x, y, boundary=boundary, normalize_kernel=normalize_kernel) 353 354 # boundary, normalize_kernel == False 355 rslt = { 356 (None): [0., 0., 0., 0., 9., 0., 0., 0., 0.], 357 ('fill'): [-1., 3., 6., 8., 9., 10., 21., 33., 46.], 358 ('wrap'): [-31., -21., -11., -1., 9., 10., 20., 30., 40.], 359 ('extend'): [-5., 0., 4., 7., 9., 10., 12., 15., 19.] 360 }[boundary] 361 362 assert_array_almost_equal_nulp(z, np.array(rslt, dtype='>f8'), 10) 363 364 @pytest.mark.parametrize('preserve_nan', PRESERVE_NAN_OPTIONS) 365 def test_int_masked_array(self, preserve_nan): 366 """ 367 Test that convolve works correctly with integer masked arrays. 368 """ 369 370 x = ma.array([3, 5, 7, 11, 13], mask=[0, 0, 1, 0, 0], fill_value=0.) 371 y = np.array([1., 1., 1.], dtype='>f8') 372 373 z = convolve(x, y, preserve_nan=preserve_nan) 374 375 if preserve_nan: 376 assert np.isnan(z[2]) 377 z[2] = 8 378 379 assert_array_almost_equal_nulp(z, (8/3., 4, 8, 12, 8), 10) 380 381 382class TestConvolve2D: 383 def test_list(self): 384 """ 385 Test that convolve works correctly when inputs are lists 386 """ 387 x = [[1, 1, 1], 388 [1, 1, 1], 389 [1, 1, 1]] 390 391 z = convolve(x, x, boundary='fill', fill_value=1, normalize_kernel=True) 392 assert_array_almost_equal_nulp(z, x, 10) 393 z = convolve(x, x, boundary='fill', fill_value=1, normalize_kernel=False) 394 assert_array_almost_equal_nulp(z, np.array(x, float)*9, 10) 395 396 @pytest.mark.parametrize(('dtype_array', 'dtype_kernel'), VALID_DTYPE_MATRIX) 397 def test_dtype(self, dtype_array, dtype_kernel): 398 ''' 399 Test that 32- and 64-bit floats are correctly handled 400 ''' 401 402 x = np.array([[1., 2., 3.], 403 [4., 5., 6.], 404 [7., 8., 9.]], dtype=dtype_array) 405 406 y = np.array([[0., 0., 0.], 407 [0., 1., 0.], 408 [0., 0., 0.]], dtype=dtype_kernel) 409 410 z = convolve(x, y) 411 412 assert x.dtype == z.dtype 413 414 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 415 def test_unity_1x1_none(self, boundary): 416 ''' 417 Test that a 1x1 unit kernel returns the same array 418 ''' 419 420 x = np.array([[1., 2., 3.], 421 [4., 5., 6.], 422 [7., 8., 9.]], dtype='>f8') 423 424 y = np.array([[1.]], dtype='>f8') 425 426 z = convolve(x, y, boundary=boundary) 427 428 assert np.all(z == x) 429 430 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 431 def test_unity_3x3(self, boundary): 432 ''' 433 Test that a 3x3 unit kernel returns the same array (except when 434 boundary is None). 435 ''' 436 437 x = np.array([[1., 2., 3.], 438 [4., 5., 6.], 439 [7., 8., 9.]], dtype='>f8') 440 441 y = np.array([[0., 0., 0.], 442 [0., 1., 0.], 443 [0., 0., 0.]], dtype='>f8') 444 445 z = convolve(x, y, boundary=boundary) 446 447 if boundary is None: 448 assert np.all(z == np.array([[0., 0., 0.], 449 [0., 5., 0.], 450 [0., 0., 0.]], dtype='>f8')) 451 else: 452 assert np.all(z == x) 453 454 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 455 def test_uniform_3x3(self, boundary): 456 ''' 457 Test that the different modes are producing the correct results using 458 a 3x3 uniform kernel. 459 ''' 460 461 x = np.array([[0., 0., 3.], 462 [1., 0., 0.], 463 [0., 2., 0.]], dtype='>f8') 464 465 y = np.array([[1., 1., 1.], 466 [1., 1., 1.], 467 [1., 1., 1.]], dtype='>f8') 468 469 z = convolve(x, y, boundary=boundary, normalize_kernel=False) 470 471 if boundary is None: 472 assert_array_almost_equal_nulp(z, np.array([[0., 0., 0.], 473 [0., 6., 0.], 474 [0., 0., 0.]], dtype='>f8'), 10) 475 elif boundary == 'fill': 476 assert_array_almost_equal_nulp(z, np.array([[1., 4., 3.], 477 [3., 6., 5.], 478 [3., 3., 2.]], dtype='>f8'), 10) 479 elif boundary == 'wrap': 480 assert_array_almost_equal_nulp(z, np.array([[6., 6., 6.], 481 [6., 6., 6.], 482 [6., 6., 6.]], dtype='>f8'), 10) 483 else: 484 assert_array_almost_equal_nulp(z, np.array([[2., 7., 12.], 485 [4., 6., 8.], 486 [6., 5., 4.]], dtype='>f8'), 10) 487 488 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 489 def test_unity_3x3_withnan(self, boundary): 490 ''' 491 Test that a 3x3 unit kernel returns the same array (except when 492 boundary is None). This version includes a NaN value in the original 493 array. 494 ''' 495 496 x = np.array([[1., 2., 3.], 497 [4., np.nan, 6.], 498 [7., 8., 9.]], dtype='>f8') 499 500 y = np.array([[0., 0., 0.], 501 [0., 1., 0.], 502 [0., 0., 0.]], dtype='>f8') 503 504 z = convolve(x, y, boundary=boundary, nan_treatment='fill', 505 preserve_nan=True) 506 507 assert np.isnan(z[1, 1]) 508 x = np.nan_to_num(z) 509 z = np.nan_to_num(z) 510 511 if boundary is None: 512 assert np.all(z == np.array([[0., 0., 0.], 513 [0., 0., 0.], 514 [0., 0., 0.]], dtype='>f8')) 515 else: 516 assert np.all(z == x) 517 518 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 519 def test_uniform_3x3_withnanfilled(self, boundary): 520 ''' 521 Test that the different modes are producing the correct results using 522 a 3x3 uniform kernel. This version includes a NaN value in the 523 original array. 524 ''' 525 526 x = np.array([[0., 0., 4.], 527 [1., np.nan, 0.], 528 [0., 3., 0.]], dtype='>f8') 529 530 y = np.array([[1., 1., 1.], 531 [1., 1., 1.], 532 [1., 1., 1.]], dtype='>f8') 533 534 z = convolve(x, y, boundary=boundary, nan_treatment='fill', 535 normalize_kernel=False) 536 537 if boundary is None: 538 assert_array_almost_equal_nulp(z, np.array([[0., 0., 0.], 539 [0., 8., 0.], 540 [0., 0., 0.]], dtype='>f8'), 10) 541 elif boundary == 'fill': 542 assert_array_almost_equal_nulp(z, np.array([[1., 5., 4.], 543 [4., 8., 7.], 544 [4., 4., 3.]], dtype='>f8'), 10) 545 elif boundary == 'wrap': 546 assert_array_almost_equal_nulp(z, np.array([[8., 8., 8.], 547 [8., 8., 8.], 548 [8., 8., 8.]], dtype='>f8'), 10) 549 elif boundary == 'extend': 550 assert_array_almost_equal_nulp(z, np.array([[2., 9., 16.], 551 [5., 8., 11.], 552 [8., 7., 6.]], dtype='>f8'), 10) 553 else: 554 raise ValueError("Invalid boundary specification") 555 556 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 557 def test_uniform_3x3_withnaninterped(self, boundary): 558 ''' 559 Test that the different modes are producing the correct results using 560 a 3x3 uniform kernel. This version includes a NaN value in the 561 original array. 562 ''' 563 564 x = np.array([[0., 0., 4.], 565 [1., np.nan, 0.], 566 [0., 3., 0.]], dtype='>f8') 567 568 y = np.array([[1., 1., 1.], 569 [1., 1., 1.], 570 [1., 1., 1.]], dtype='>f8') 571 572 z = convolve(x, y, boundary=boundary, nan_treatment='interpolate', 573 normalize_kernel=True) 574 575 if boundary is None: 576 assert_array_almost_equal_nulp(z, np.array([[0., 0., 0.], 577 [0., 1., 0.], 578 [0., 0., 0.]], dtype='>f8'), 10) 579 elif boundary == 'fill': 580 assert_array_almost_equal_nulp(z, np.array([[1./8, 5./8, 4./8], 581 [4./8, 8./8, 7./8], 582 [4./8, 4./8, 3./8]], dtype='>f8'), 10) 583 elif boundary == 'wrap': 584 assert_array_almost_equal_nulp(z, np.array([[1., 1., 1.], 585 [1., 1., 1.], 586 [1., 1., 1.]], dtype='>f8'), 10) 587 elif boundary == 'extend': 588 assert_array_almost_equal_nulp(z, np.array([[2./8, 9./8, 16./8], 589 [5./8, 8./8, 11./8], 590 [8./8, 7./8, 6./8]], dtype='>f8'), 10) 591 else: 592 raise ValueError("Invalid boundary specification") 593 594 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 595 def test_non_normalized_kernel_2D(self, boundary): 596 597 x = np.array([[0., 0., 4.], 598 [1., 2., 0.], 599 [0., 3., 0.]], dtype='float') 600 601 y = np.array([[1., -1., 1.], 602 [-1., 0., -1.], 603 [1., -1., 1.]], dtype='float') 604 605 z = convolve(x, y, boundary=boundary, nan_treatment='fill', 606 normalize_kernel=False) 607 608 if boundary is None: 609 assert_array_almost_equal_nulp(z, np.array([[0., 0., 0.], 610 [0., 0., 0.], 611 [0., 0., 0.]], dtype='float'), 10) 612 elif boundary == 'fill': 613 assert_array_almost_equal_nulp(z, np.array([[1., -5., 2.], 614 [1., 0., -3.], 615 [-2., -1., -1.]], dtype='float'), 10) 616 elif boundary == 'wrap': 617 assert_array_almost_equal_nulp(z, np.array([[0., -8., 6.], 618 [5., 0., -4.], 619 [2., 3., -4.]], dtype='float'), 10) 620 elif boundary == 'extend': 621 assert_array_almost_equal_nulp(z, np.array([[2., -1., -2.], 622 [0., 0., 1.], 623 [2., -4., 2.]], dtype='float'), 10) 624 else: 625 raise ValueError("Invalid boundary specification") 626 627 628class TestConvolve3D: 629 def test_list(self): 630 """ 631 Test that convolve works correctly when inputs are lists 632 """ 633 x = [[[1, 1, 1], 634 [1, 1, 1], 635 [1, 1, 1]], 636 [[1, 1, 1], 637 [1, 1, 1], 638 [1, 1, 1]], 639 [[1, 1, 1], 640 [1, 1, 1], 641 [1, 1, 1]]] 642 643 z = convolve(x, x, boundary='fill', fill_value=1, normalize_kernel=False) 644 assert_array_almost_equal_nulp(z / 27, x, 10) 645 646 @pytest.mark.parametrize(('dtype_array', 'dtype_kernel'), VALID_DTYPE_MATRIX) 647 def test_dtype(self, dtype_array, dtype_kernel): 648 ''' 649 Test that 32- and 64-bit floats are correctly handled 650 ''' 651 652 x = np.array([[1., 2., 3.], 653 [4., 5., 6.], 654 [7., 8., 9.]], dtype=dtype_array) 655 656 y = np.array([[0., 0., 0.], 657 [0., 1., 0.], 658 [0., 0., 0.]], dtype=dtype_kernel) 659 660 z = convolve(x, y) 661 662 assert x.dtype == z.dtype 663 664 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 665 def test_unity_1x1x1_none(self, boundary): 666 ''' 667 Test that a 1x1x1 unit kernel returns the same array 668 ''' 669 670 x = np.array([[[1., 2., 1.], [2., 3., 1.], [3., 2., 5.]], 671 [[4., 3., 1.], [5., 0., 2.], [6., 1., 1.]], 672 [[7., 0., 2.], [8., 2., 3.], [9., 2., 2.]]], dtype='>f8') 673 674 y = np.array([[[1.]]], dtype='>f8') 675 676 z = convolve(x, y, boundary=boundary) 677 678 assert np.all(z == x) 679 680 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 681 def test_unity_3x3x3(self, boundary): 682 ''' 683 Test that a 3x3x3 unit kernel returns the same array (except when 684 boundary is None). 685 ''' 686 687 x = np.array([[[1., 2., 1.], [2., 3., 1.], [3., 2., 5.]], 688 [[4., 3., 1.], [5., 3., 2.], [6., 1., 1.]], 689 [[7., 0., 2.], [8., 2., 3.], [9., 2., 2.]]], dtype='>f8') 690 691 y = np.zeros((3, 3, 3), dtype='>f8') 692 y[1, 1, 1] = 1. 693 694 z = convolve(x, y, boundary=boundary) 695 696 if boundary is None: 697 assert np.all(z == np.array([[[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], 698 [[0., 0., 0.], [0., 3., 0.], [0., 0., 0.]], 699 [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]], dtype='>f8')) 700 else: 701 assert np.all(z == x) 702 703 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 704 def test_uniform_3x3x3(self, boundary): 705 ''' 706 Test that the different modes are producing the correct results using 707 a 3x3 uniform kernel. 708 ''' 709 710 x = np.array([[[1., 2., 1.], [2., 3., 1.], [3., 2., 5.]], 711 [[4., 3., 1.], [5., 3., 2.], [6., 1., 1.]], 712 [[7., 0., 2.], [8., 2., 3.], [9., 2., 2.]]], dtype='>f8') 713 714 y = np.ones((3, 3, 3), dtype='>f8') 715 716 z = convolve(x, y, boundary=boundary, normalize_kernel=False) 717 718 if boundary is None: 719 assert_array_almost_equal_nulp(z, np.array([[[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], 720 [[0., 0., 0.], [0., 81., 0.], [0., 0., 0.]], 721 [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]], dtype='>f8'), 10) 722 elif boundary == 'fill': 723 assert_array_almost_equal_nulp(z, np.array([[[23., 28., 16.], [35., 46., 25.], [25., 34., 18.]], 724 [[40., 50., 23.], [63., 81., 36.], [46., 60., 27.]], 725 [[32., 40., 16.], [50., 61., 22.], [36., 44., 16.]]], dtype='>f8'), 10) 726 elif boundary == 'wrap': 727 assert_array_almost_equal_nulp(z, np.array([[[81., 81., 81.], [81., 81., 81.], [81., 81., 81.]], 728 [[81., 81., 81.], [81., 81., 81.], [81., 81., 81.]], 729 [[81., 81., 81.], [81., 81., 81.], [81., 81., 81.]]], dtype='>f8'), 10) 730 else: 731 assert_array_almost_equal_nulp(z, np.array([[[65., 54., 43.], [75., 66., 57.], [85., 78., 71.]], 732 [[96., 71., 46.], [108., 81., 54.], [120., 91., 62.]], 733 [[127., 88., 49.], [141., 96., 51.], [155., 104., 53.]]], dtype='>f8'), 10) 734 735 @pytest.mark.parametrize(('boundary', 'nan_treatment'), 736 itertools.product(BOUNDARY_OPTIONS, 737 NANHANDLING_OPTIONS)) 738 def test_unity_3x3x3_withnan(self, boundary, nan_treatment): 739 ''' 740 Test that a 3x3x3 unit kernel returns the same array (except when 741 boundary is None). This version includes a NaN value in the original 742 array. 743 ''' 744 745 x = np.array([[[1., 2., 1.], [2., 3., 1.], [3., 2., 5.]], 746 [[4., 3., 1.], [5., np.nan, 2.], [6., 1., 1.]], 747 [[7., 0., 2.], [8., 2., 3.], [9., 2., 2.]]], dtype='>f8') 748 749 y = np.zeros((3, 3, 3), dtype='>f8') 750 y[1, 1, 1] = 1. 751 752 z = convolve(x, y, boundary=boundary, nan_treatment=nan_treatment, 753 preserve_nan=True) 754 755 assert np.isnan(z[1, 1, 1]) 756 x = np.nan_to_num(z) 757 z = np.nan_to_num(z) 758 759 if boundary is None: 760 assert np.all(z == np.array([[[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], 761 [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], 762 [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]], dtype='>f8')) 763 else: 764 assert np.all(z == x) 765 766 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 767 def test_uniform_3x3x3_withnan_filled(self, boundary): 768 ''' 769 Test that the different modes are producing the correct results using 770 a 3x3 uniform kernel. This version includes a NaN value in the 771 original array. 772 ''' 773 774 x = np.array([[[1., 2., 1.], [2., 3., 1.], [3., 2., 5.]], 775 [[4., 3., 1.], [5., np.nan, 2.], [6., 1., 1.]], 776 [[7., 0., 2.], [8., 2., 3.], [9., 2., 2.]]], dtype='>f8') 777 778 y = np.ones((3, 3, 3), dtype='>f8') 779 780 z = convolve(x, y, boundary=boundary, nan_treatment='fill', 781 normalize_kernel=False) 782 783 if boundary is None: 784 assert_array_almost_equal_nulp(z, np.array([[[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], 785 [[0., 0., 0.], [0., 78., 0.], [0., 0., 0.]], 786 [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]], dtype='>f8'), 10) 787 elif boundary == 'fill': 788 assert_array_almost_equal_nulp(z, np.array([[[20., 25., 13.], 789 [32., 43., 22.], 790 [22., 31., 15.]], 791 [[37., 47., 20.], 792 [60., 78., 33.], 793 [43., 57., 24.]], 794 [[29., 37., 13.], 795 [47., 58., 19.], 796 [33., 41., 13.]]], dtype='>f8'), 10) 797 elif boundary == 'wrap': 798 assert_array_almost_equal_nulp(z, np.array([[[78., 78., 78.], [78., 78., 78.], [78., 78., 78.]], 799 [[78., 78., 78.], [78., 78., 78.], [78., 78., 78.]], 800 [[78., 78., 78.], [78., 78., 78.], [78., 78., 78.]]], dtype='>f8'), 10) 801 elif boundary == 'extend': 802 assert_array_almost_equal_nulp(z, np.array([[[62., 51., 40.], 803 [72., 63., 54.], 804 [82., 75., 68.]], 805 [[93., 68., 43.], 806 [105., 78., 51.], 807 [117., 88., 59.]], 808 [[124., 85., 46.], 809 [138., 93., 48.], 810 [152., 101., 50.]]], 811 dtype='>f8'), 10) 812 else: 813 raise ValueError("Invalid Boundary Option") 814 815 @pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 816 def test_uniform_3x3x3_withnan_interped(self, boundary): 817 ''' 818 Test that the different modes are producing the correct results using 819 a 3x3 uniform kernel. This version includes a NaN value in the 820 original array. 821 ''' 822 823 x = np.array([[[1., 2., 1.], [2., 3., 1.], [3., 2., 5.]], 824 [[4., 3., 1.], [5., np.nan, 2.], [6., 1., 1.]], 825 [[7., 0., 2.], [8., 2., 3.], [9., 2., 2.]]], dtype='>f8') 826 827 y = np.ones((3, 3, 3), dtype='>f8') 828 829 z = convolve(x, y, boundary=boundary, nan_treatment='interpolate', 830 normalize_kernel=True) 831 832 kernsum = y.sum() - 1 # one nan is missing 833 mid = x[np.isfinite(x)].sum() / kernsum 834 835 if boundary is None: 836 assert_array_almost_equal_nulp(z, np.array([[[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], 837 [[0., 0., 0.], [0., 78., 0.], [0., 0., 0.]], 838 [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]], 839 dtype='>f8')/kernsum, 10) 840 elif boundary == 'fill': 841 assert_array_almost_equal_nulp(z, np.array([[[20., 25., 13.], 842 [32., 43., 22.], 843 [22., 31., 15.]], 844 [[37., 47., 20.], 845 [60., 78., 33.], 846 [43., 57., 24.]], 847 [[29., 37., 13.], 848 [47., 58., 19.], 849 [33., 41., 13.]]], 850 dtype='>f8')/kernsum, 10) 851 elif boundary == 'wrap': 852 assert_array_almost_equal_nulp(z, np.tile(mid.astype('>f8'), [3, 3, 3]), 10) 853 elif boundary == 'extend': 854 assert_array_almost_equal_nulp(z, np.array([[[62., 51., 40.], 855 [72., 63., 54.], 856 [82., 75., 68.]], 857 [[93., 68., 43.], 858 [105., 78., 51.], 859 [117., 88., 59.]], 860 [[124., 85., 46.], 861 [138., 93., 48.], 862 [152., 101., 50.]]], 863 dtype='>f8')/kernsum, 10) 864 else: 865 raise ValueError("Invalid Boundary Option") 866 867 868@pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 869def test_asymmetric_kernel(boundary): 870 ''' 871 Regression test for #6264: make sure that asymmetric convolution 872 functions go the right direction 873 ''' 874 875 x = np.array([3., 0., 1.], dtype='>f8') 876 877 y = np.array([1, 2, 3], dtype='>f8') 878 879 z = convolve(x, y, boundary=boundary, normalize_kernel=False) 880 881 if boundary == 'fill': 882 assert_array_almost_equal_nulp(z, np.array([6., 10., 2.], dtype='float'), 10) 883 elif boundary is None: 884 assert_array_almost_equal_nulp(z, np.array([0., 10., 0.], dtype='float'), 10) 885 elif boundary == 'extend': 886 assert_array_almost_equal_nulp(z, np.array([15., 10., 3.], dtype='float'), 10) 887 elif boundary == 'wrap': 888 assert_array_almost_equal_nulp(z, np.array([9., 10., 5.], dtype='float'), 10) 889 890 891@pytest.mark.parametrize('ndims', (1, 2, 3)) 892def test_convolution_consistency(ndims): 893 894 np.random.seed(0) 895 array = np.random.randn(*([3]*ndims)) 896 np.random.seed(0) 897 kernel = np.random.rand(*([3]*ndims)) 898 899 conv_f = convolve_fft(array, kernel, boundary='fill') 900 conv_d = convolve(array, kernel, boundary='fill') 901 902 assert_array_almost_equal_nulp(conv_f, conv_d, 30) 903 904 905def test_astropy_convolution_against_numpy(): 906 x = np.array([1, 2, 3]) 907 y = np.array([5, 4, 3, 2, 1]) 908 909 assert_array_almost_equal(np.convolve(y, x, 'same'), 910 convolve(y, x, normalize_kernel=False)) 911 assert_array_almost_equal(np.convolve(y, x, 'same'), 912 convolve_fft(y, x, normalize_kernel=False)) 913 914 915@pytest.mark.skipif('not HAS_SCIPY') 916def test_astropy_convolution_against_scipy(): 917 from scipy.signal import fftconvolve 918 x = np.array([1, 2, 3]) 919 y = np.array([5, 4, 3, 2, 1]) 920 921 assert_array_almost_equal(fftconvolve(y, x, 'same'), 922 convolve(y, x, normalize_kernel=False)) 923 assert_array_almost_equal(fftconvolve(y, x, 'same'), 924 convolve_fft(y, x, normalize_kernel=False)) 925 926 927@pytest.mark.skipif('not HAS_PANDAS') 928def test_regression_6099(): 929 import pandas 930 931 wave = np.array(np.linspace(5000, 5100, 10)) 932 boxcar = 3 933 nonseries_result = convolve(wave, np.ones((boxcar,))/boxcar) 934 935 wave_series = pandas.Series(wave) 936 series_result = convolve(wave_series, np.ones((boxcar,))/boxcar) 937 938 assert_array_almost_equal(nonseries_result, series_result) 939 940 941def test_invalid_array_convolve(): 942 kernel = np.ones(3)/3. 943 944 with pytest.raises(TypeError): 945 convolve('glork', kernel) 946 947 948@pytest.mark.parametrize(('boundary'), BOUNDARY_OPTIONS) 949def test_non_square_kernel_asymmetric(boundary): 950 # Regression test for a bug that occurred when using non-square kernels in 951 # 2D when using boundary=None 952 kernel = np.array([[1, 2, 3, 2, 1], [0, 1, 2, 1, 0], [0, 0, 0, 0, 0]]) 953 image = np.zeros((13, 13)) 954 image[6, 6] = 1 955 result = convolve(image, kernel, normalize_kernel=False, boundary=boundary) 956 assert_allclose(result[5:8, 4:9], kernel) 957 958 959@pytest.mark.parametrize(('boundary', 'normalize_kernel'), 960 itertools.product(BOUNDARY_OPTIONS, 961 NORMALIZE_OPTIONS)) 962def test_uninterpolated_nan_regions(boundary, normalize_kernel): 963 #8086 964 # Test NaN interpolation of contiguous NaN regions with kernels of size 965 # identical and greater than that of the region of NaN values. 966 967 # Test case: kernel.shape == NaN_region.shape 968 kernel = Gaussian2DKernel(1, 5, 5) 969 nan_centroid = np.full(kernel.shape, np.nan) 970 image = np.pad(nan_centroid, pad_width=kernel.shape[0]*2, mode='constant', 971 constant_values=1) 972 with pytest.warns(AstropyUserWarning, 973 match=r"nan_treatment='interpolate', however, NaN values detected " 974 r"post convolution. A contiguous region of NaN values, larger " 975 r"than the kernel size, are present in the input array. " 976 r"Increase the kernel size to avoid this."): 977 result = convolve(image, kernel, boundary=boundary, nan_treatment='interpolate', 978 normalize_kernel=normalize_kernel) 979 assert(np.any(np.isnan(result))) 980 981 # Test case: kernel.shape > NaN_region.shape 982 nan_centroid = np.full((kernel.shape[0]-1, kernel.shape[1]-1), np.nan) # 1 smaller than kerenel 983 image = np.pad(nan_centroid, pad_width=kernel.shape[0]*2, mode='constant', 984 constant_values=1) 985 result = convolve(image, kernel, boundary=boundary, nan_treatment='interpolate', 986 normalize_kernel=normalize_kernel) 987 assert(~np.any(np.isnan(result))) # Note: negation 988 989 990def test_regressiontest_issue9168(): 991 """ 992 Issue #9168 pointed out that kernels can be (unitless) quantities, which 993 leads to crashes when inplace modifications are made to arrays in 994 convolve/convolve_fft, so we now strip the quantity aspects off of kernels. 995 """ 996 997 x = np.array([[1., 2., 3.], 998 [4., 5., 6.], 999 [7., 8., 9.]],) 1000 1001 kernel_fwhm = 1*u.arcsec 1002 pixel_size = 1*u.arcsec 1003 1004 kernel = Gaussian2DKernel(x_stddev=kernel_fwhm/pixel_size) 1005 1006 result = convolve_fft(x, kernel, boundary='fill', fill_value=np.nan, 1007 preserve_nan=True) 1008 result = convolve(x, kernel, boundary='fill', fill_value=np.nan, 1009 preserve_nan=True) 1010 1011 1012def test_convolve_nan_zero_sum_kernel(): 1013 with pytest.raises(ValueError, 1014 match="Setting nan_treatment='interpolate' " 1015 "requires the kernel to be normalized, but the " 1016 "input kernel has a sum close to zero. For a " 1017 "zero-sum kernel and data with NaNs, set " 1018 "nan_treatment='fill'."): 1019 convolve([1, np.nan, 3], [-1, 2, -1], normalize_kernel=False) 1020