# Licensed under a 3-clause BSD style license - see LICENSE.rst

import itertools
from contextlib import nullcontext

import pytest
import numpy as np
from numpy.testing import assert_allclose, assert_array_equal, assert_array_almost_equal_nulp

from astropy.convolution.convolve import convolve_fft, convolve
from astropy.utils.exceptions import AstropyUserWarning
from astropy import units as u

VALID_DTYPES = ('>f4', '<f4', '>f8', '<f8')
VALID_DTYPE_MATRIX = list(itertools.product(VALID_DTYPES, VALID_DTYPES))

BOUNDARY_OPTIONS = [None, 'fill', 'wrap']
NANTREATMENT_OPTIONS = ('interpolate', 'fill')
NORMALIZE_OPTIONS = [True, False]
PRESERVE_NAN_OPTIONS = [True, False]

"""
What does convolution mean?  We use the 'same size' assumption here (i.e.,
you expect an array of the exact same size as the one you put in)
Convolving any array with a kernel that is [1] should result in the same array returned
Working example array: [1, 2, 3, 4, 5]
Convolved with [1] = [1, 2, 3, 4, 5]
Convolved with [1, 1] = [1, 3, 5, 7, 9] THIS IS NOT CONSISTENT!
Convolved with [1, 0] = [1, 2, 3, 4, 5]
Convolved with [0, 1] = [0, 1, 2, 3, 4]
"""

# NOTE: use_numpy_fft is redundant if you don't have FFTW installed
option_names = ('boundary', 'nan_treatment', 'normalize_kernel')
options = list(itertools.product(BOUNDARY_OPTIONS,
                                 NANTREATMENT_OPTIONS,
                                 (True, False),
                                 ))
option_names_preserve_nan = ('boundary', 'nan_treatment',
                             'normalize_kernel', 'preserve_nan')
options_preserve_nan = list(itertools.product(BOUNDARY_OPTIONS,
                                              NANTREATMENT_OPTIONS,
                                              (True, False),
                                              (True, False)))


def expected_boundary_warning(boundary=None):
    # Helper that returns the appropriate context manager for the boundary=None
    # warning depending on the value of boundary.
    if boundary is None:
        ctx = pytest.warns(AstropyUserWarning,
                           match='The convolve_fft version of boundary=None '
                                 'is equivalent to the convolve boundary=\'fill\'')
    else:
        ctx = nullcontext()
    return ctx


def assert_floatclose(x, y):
    """Assert arrays are close to within expected floating point rounding.

    Check that the result is correct at the precision expected for 64 bit
    numbers, taking account that the tolerance has to reflect that all powers
    in the FFTs enter our values.
    """
    # The number used is set by the fact that the Windows FFT sometimes
    # returns an answer that is EXACTLY 10*np.spacing.
    assert_allclose(x, y, atol=10*np.spacing(x.max()), rtol=0.)


class TestConvolve1D:

    @pytest.mark.parametrize(option_names, options)
    def test_quantity(self, boundary, nan_treatment, normalize_kernel):
        """
        Test that convolve_fft works correctly when input array is a Quantity
        """

        x = np.array([1., 4., 5., 6., 5., 7., 8.], dtype='float64') * u.ph
        y = np.array([0.2, 0.6, 0.2], dtype='float64')

        with expected_boundary_warning(boundary=boundary):
            z = convolve_fft(x, y, boundary=boundary,
                             nan_treatment=nan_treatment,
                             normalize_kernel=normalize_kernel)

        assert x.unit == z.unit

    @pytest.mark.parametrize(option_names, options)
    def test_unity_1_none(self, boundary, nan_treatment, normalize_kernel):
        '''
        Test that a unit kernel with a single element returns the same array
        '''

        x = np.array([1., 2., 3.], dtype='float64')

        y = np.array([1.], dtype='float64')

        with expected_boundary_warning(boundary=boundary):
            z = convolve_fft(x, y, boundary=boundary,
                             nan_treatment=nan_treatment,
                             normalize_kernel=normalize_kernel)

        assert_floatclose(z, x)

    @pytest.mark.parametrize(option_names, options)
    def test_unity_3(self, boundary, nan_treatment, normalize_kernel):
        '''
        Test that a unit kernel with three elements returns the same array
        (except when boundary is None).
        '''

        x = np.array([1., 2., 3.], dtype='float64')

        y = np.array([0., 1., 0.], dtype='float64')

        with expected_boundary_warning(boundary=boundary):
            z = convolve_fft(x, y, boundary=boundary,
                             nan_treatment=nan_treatment,
                             normalize_kernel=normalize_kernel)

        assert_floatclose(z, x)

    @pytest.mark.parametrize(option_names, options)
    def test_uniform_3(self, boundary, nan_treatment, normalize_kernel):
        '''
        Test that the different modes are producing the correct results using
        a uniform kernel with three elements
        '''

        x = np.array([1., 0., 3.], dtype='float64')

        y = np.array([1., 1., 1.], dtype='float64')

        with expected_boundary_warning(boundary=boundary):
            z = convolve_fft(x, y, boundary=boundary,
                             nan_treatment=nan_treatment,
                             normalize_kernel=normalize_kernel)

        answer_key = (boundary, nan_treatment, normalize_kernel)

        answer_dict = {
            'sum_fill_zeros': np.array([1., 4., 3.], dtype='float64'),
            'average_fill_zeros': np.array([1 / 3., 4 / 3., 1.], dtype='float64'),
            'sum_wrap': np.array([4., 4., 4.], dtype='float64'),
            'average_wrap': np.array([4 / 3., 4 / 3., 4 / 3.], dtype='float64'),
        }

        result_dict = {
            # boundary, nan_treatment, normalize_kernel
            ('fill', 'interpolate', True): answer_dict['average_fill_zeros'],
            ('wrap', 'interpolate', True): answer_dict['average_wrap'],
            ('fill', 'interpolate', False): answer_dict['sum_fill_zeros'],
            ('wrap', 'interpolate', False): answer_dict['sum_wrap'],
        }
        for k in list(result_dict.keys()):
            result_dict[(k[0], 'fill', k[2])] = result_dict[k]
        for k in list(result_dict.keys()):
            if k[0] == 'fill':
                result_dict[(None, k[1], k[2])] = result_dict[k]

        assert_floatclose(z, result_dict[answer_key])

    @pytest.mark.parametrize(option_names, options)
    def test_halfity_3(self, boundary, nan_treatment, normalize_kernel):
        '''
        Test that the different modes are producing the correct results using
        a uniform, non-unity kernel with three elements
        '''

        x = np.array([1., 0., 3.], dtype='float64')

        y = np.array([0.5, 0.5, 0.5], dtype='float64')

        with expected_boundary_warning(boundary=boundary):
            z = convolve_fft(x, y, boundary=boundary,
                             nan_treatment=nan_treatment,
                             normalize_kernel=normalize_kernel)

        answer_dict = {
            'sum': np.array([0.5, 2.0, 1.5], dtype='float64'),
            'sum_zeros': np.array([0.5, 2., 1.5], dtype='float64'),
            'sum_nozeros': np.array([0.5, 2., 1.5], dtype='float64'),
            'average': np.array([1 / 3., 4 / 3., 1.], dtype='float64'),
            'sum_wrap': np.array([2., 2., 2.], dtype='float64'),
            'average_wrap': np.array([4 / 3., 4 / 3., 4 / 3.], dtype='float64'),
            'average_zeros': np.array([1 / 3., 4 / 3., 1.], dtype='float64'),
            'average_nozeros': np.array([0.5, 4 / 3., 1.5], dtype='float64'),
        }

        if normalize_kernel:
            answer_key = 'average'
        else:
            answer_key = 'sum'

        if boundary == 'wrap':
            answer_key += '_wrap'
        else:
            # average = average_zeros; sum = sum_zeros
            answer_key += '_zeros'

        assert_floatclose(z, answer_dict[answer_key])

    @pytest.mark.parametrize(option_names_preserve_nan, options_preserve_nan)
    def test_unity_3_withnan(self, boundary, nan_treatment, normalize_kernel,
                             preserve_nan):
        '''
        Test that a unit kernel with three elements returns the same array
        (except when boundary is None). This version includes a NaN value in
        the original array.
        '''

        x = np.array([1., np.nan, 3.], dtype='float64')

        y = np.array([0., 1., 0.], dtype='float64')

        with expected_boundary_warning(boundary=boundary):
            z = convolve_fft(x, y, boundary=boundary,
                             nan_treatment=nan_treatment,
                             normalize_kernel=normalize_kernel,
                             preserve_nan=preserve_nan)

        if preserve_nan:
            assert np.isnan(z[1])

        z = np.nan_to_num(z)

        assert_floatclose(z, [1., 0., 3.])

    inputs = (np.array([1., np.nan, 3.], dtype='float64'),
              np.array([1., np.inf, 3.], dtype='float64'))
    outputs = (np.array([1., 0., 3.], dtype='float64'),
               np.array([1., 0., 3.], dtype='float64'))
    options_unity1withnan = list(itertools.product(BOUNDARY_OPTIONS,
                                                   NANTREATMENT_OPTIONS,
                                                   (True, False),
                                                   (True, False),
                                                   inputs, outputs))

    @pytest.mark.parametrize(option_names_preserve_nan + ('inval', 'outval'),
                             options_unity1withnan)
    def test_unity_1_withnan(self, boundary, nan_treatment, normalize_kernel,
                             preserve_nan, inval, outval):
        '''
        Test that a unit kernel with three elements returns the same array
        (except when boundary is None). This version includes a NaN value in
        the original array.
        '''

        x = inval

        y = np.array([1.], dtype='float64')

        with expected_boundary_warning(boundary=boundary):
            z = convolve_fft(x, y, boundary=boundary,
                             nan_treatment=nan_treatment,
                             normalize_kernel=normalize_kernel,
                             preserve_nan=preserve_nan)

        if preserve_nan:
            assert np.isnan(z[1])

        z = np.nan_to_num(z)

        assert_floatclose(z, outval)

    @pytest.mark.parametrize(option_names_preserve_nan, options_preserve_nan)
    def test_uniform_3_withnan(self, boundary, nan_treatment,
                               normalize_kernel, preserve_nan):
        '''
        Test that the different modes are producing the correct results using