[Bug] Testing issue: py-elephant (JUSUF, Galileo100)

[Moritz-Alexander-Kern]

Describe the bug
The ebrains TC team have notified me about the following issue when testing Elephant on JUSUF and Galileo100.

When trying to install elephant on JUSUF, two tests fail: (on Galileo100 only the WLPI ground truth test fails)

elephant/test/test_phase_analysis.py::WeightedPhaseLagIndexTestCase::test_WPLI_ground_truth_consistency_real_LFP_dataset
elephant/test/test_spectral.py::MultitaperCoherenceTestCase::test_multitaper_cohere_perfect_cohere

To Reproduce

1. $ spack spack install --test root py-elephant

Expected behavior
No error should be raised when running the unit tests.

Environment
$ spack debug report

• Spack: 0.18.1 (57226fd1f4c6bc8c55fe0bfde9450a326e957f0d)
• Python: 3.9.12
• Platform: linux-centos8-cascadelake
• Concretizer: clingo
• tqdm-4.62.3
• statsmodels-0.13.1
• scikit-learn-1.1.1
• pytest-6.2.5
• pandas-1.4.2
• six-1.16.0
• neo-0.12.0
• scipy-1.8.1
• quantities-0.14.1
• numpy-1.21.6
• jinja2-3.0.3

Error

=================================== FAILURES ===================================
_ WeightedPhaseLagIndexTestCase.test_WPLI_ground_truth_consistency_real_LFP_dataset _

self = <elephant.test.test_phase_analysis.WeightedPhaseLagIndexTestCase testMethod=test_WPLI_ground_truth_consistency_real_LFP_dataset>

    def test_WPLI_ground_truth_consistency_real_LFP_dataset(self):
        """
        Test if the WPLI is consistent with the reference implementation
        ft_connectivity_wpli() of the MATLAB-package FieldTrip using
        LFP-dataset cuttings from the multielectrode-grasp  G-Node GIN
        repository, which can be found here:
        https://doi.gin.g-node.org/10.12751/g-node.f83565/
        The cutting was performed with this python-script:
        multielectrode_grasp_i140703-001_cutting_script_TS_ON_to_GO_ON.py
        which is available on https://gin.g-node.org/INM-6/elephant-data
        in folder validation/phase_analysis/weighted_phase_lag_index/scripts,
        where also the MATLAB-script for ground-truth generation is located.
        """
        # Quantity-input
        with self.subTest(msg="Quantity input"):
            freq, wpli = elephant.phase_analysis.weighted_phase_lag_index(
                self.lfps1_real, self.lfps2_real, self.sf1_real)
>           np.testing.assert_allclose(
                wpli, self.wpli_ground_truth_ft_connectivity_wpli_real,
                atol=self.tolerance, rtol=self.tolerance, equal_nan=True)
E           AssertionError: 
E           Not equal to tolerance rtol=1e-15, atol=1e-15
E           
E           Mismatched elements: 138 / 1051 (13.1%)
E           Max absolute difference: 1.33781874e-14
E           Max relative difference: 1.19137161e-11
E            x: array([     nan, 0.470965, 0.566099, ..., 0.108376, 0.069116,      nan])
E            y: array([     nan, 0.470965, 0.566099, ..., 0.108376, 0.069116,      nan])

elephant/test/test_phase_analysis.py:497: AssertionError
______ MultitaperCoherenceTestCase.test_multitaper_cohere_perfect_cohere _______

self = <elephant.test.test_spectral.MultitaperCoherenceTestCase testMethod=test_multitaper_cohere_perfect_cohere>

    def test_multitaper_cohere_perfect_cohere(self):
        # Generate dummy data
        data_length = 10000
        sampling_period = 0.001
        signal_freq = 100.0
        noise = np.random.normal(size=(1, data_length))
        time_points = np.arange(0, data_length * sampling_period,
                                sampling_period)
        signal = np.cos(2 * np.pi * signal_freq * time_points) + noise
    
        # Estimate coherence and phase lag with the multitaper method
        freq1, coh, phase_lag = elephant.spectral.multitaper_coherence(
            signal,
            signal,
            fs=1/sampling_period,
            n_segments=16)
    
        self.assertTrue((coh == np.ones(coh.size)).all())
>       self.assertTrue((phase_lag == np.zeros(phase_lag.size)).all())
E       AssertionError: False is not true

elephant/test/test_spectral.py:867: AssertionError

[Moritz-Alexander-Kern]

At first glance:

Based on the information provided, it seems that the discrepancy in the results is likely due to differences in how floating-point operations are handled on the respective hardware or compiler. These variations seem to lead to minor differences in the output.

For WLPI:
The max absolute difference of 1.33781874e-14 does not hint to a fundamental mistake in the Elephant implementation of WLPI .

It might be the tolerances set for the unit test with np.allclose are too tight, causing even small differences to be flagged as errors. Therefore the fix is likely going to be a higher tolerance for this unit test (?).

Update: Tolerance for WLPI fixed.

For MultitaperCoherence:
It's possible that this situation is similar to the one encountered in WLPI. To address it, one potential initial approach could involve modifying the assertion statement to compare the arrays using a function like np.allclose. By doing so, one would obtain an approximation of the discrepancy between the expected and calculated values in case an error is raised.

[Moritz-Alexander-Kern]

For MultitaperCoherence this might be related to: numpy/numpy#24000 .

I tracked it down to this line:

elephant/elephant/spectral.py

Lines 703 to 704 in 2651da9

	temp = (spectrum_estimates[np.newaxis, :, :, :]
	* np.conjugate(spectrum_estimates[:, np.newaxis, :, :]))

[Moritz-Alexander-Kern]

For now I suggest to wait for release of numpy 1.25.1, the issue (#24000) regarding this is scheduled for this release.

[Moritz-Alexander-Kern]

Tested with numpy==1.26.2, no changes.

[Moritz-Alexander-Kern]

Tested with numpy==1.26.4, no changes.

=================================================================================== test session starts ====================================================================================
platform linux -- Python 3.10.12, pytest-8.1.1, pluggy-1.4.0
rootdir: /home/kern/git/INM-6/elephant
plugins: anyio-4.2.0
collected 33 items                                                                                                                                                                         

elephant/test/test_spectral.py .........................F.......                                                                                                                     [100%]

========================================================================================= FAILURES =========================================================================================
____________________________________________________________ MultitaperCoherenceTestCase.test_multitaper_cohere_perfect_cohere _____________________________________________________________

self = <elephant.test.test_spectral.MultitaperCoherenceTestCase testMethod=test_multitaper_cohere_perfect_cohere>

    def test_multitaper_cohere_perfect_cohere(self):
        # Generate dummy data
        data_length = 10000
        sampling_period = 0.001
        signal_freq = 100.0
        noise = np.random.normal(size=(1, data_length))
        time_points = np.arange(0, data_length * sampling_period,
                                sampling_period)
        signal = np.cos(2 * np.pi * signal_freq * time_points) + noise
    
        # Estimate coherence and phase lag with the multitaper method
        freq1, coh, phase_lag = elephant.spectral.multitaper_coherence(
            signal,
            signal,
            fs=1/sampling_period,
            n_segments=16)
    
>       np.testing.assert_array_equal(phase_lag, np.zeros(phase_lag.size))

elephant/test/test_spectral.py:980: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

args = (<built-in function eq>, array([ 0.00000000e+00,  1.70311846e-18, -5.90941469e-18,  9.99652423e-19,
        3.00045105...., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]))
kwds = {'err_msg': '', 'header': 'Arrays are not equal', 'strict': False, 'verbose': True}

    @wraps(func)
    def inner(*args, **kwds):
        with self._recreate_cm():
>           return func(*args, **kwds)
E           AssertionError: 
E           Arrays are not equal
E           
E           Mismatched elements: 587 / 589 (99.7%)
E           Max absolute difference: 7.0909479e-18
E           Max relative difference: inf
E            x: array([ 0.000000e+00,  1.703118e-18, -5.909415e-18,  9.996524e-19,
E                   3.000451e-18, -1.319078e-19,  1.620084e-18, -1.323537e-18,
E                   7.722359e-19,  2.487218e-18,  2.764620e-18,  3.428711e-18,...
E            y: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
E                  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
E                  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,...

/usr/lib/python3.10/contextlib.py:79: AssertionError