[Feature] Add optional parallel processing support for computationally intensive functions

[feruzoripov]

Description

This PR adds opt-in parallel/concurrent execution to several computationally intensive functions in NeuroKit2. All changes are backward-compatible, default behavior remains sequential, and no new hard dependencies are introduced (joblib is already an optional dependency).

Motivation

NeuroKit2 currently runs all computations sequentially on a single core. For researchers working with long recordings, multi-channel EEG, or batch processing pipelines, this leaves significant performance on the table. Several functions contain embarrassingly parallel workloads (independent loop iterations, independent signal pipelines) that can be sped up with minimal API changes.

Proposed Changes

New parameters

Function	Parameter	What it parallelizes
bio_process()	parallel=False	Runs ECG/RSP/EDA/EMG/PPG/EOG pipelines concurrently via concurrent.futures.ProcessPoolExecutor
ecg_findpeaks() (ProMAC)	n_jobs=1	Runs 10+ peak detection methods in parallel via joblib.Parallel
entropy_multiscale()	n_jobs=1	Computes entropy at each scale factor in parallel via joblib.Parallel
eeg_power()	n_jobs=1	Processes EEG channels in parallel via joblib.Parallel
microstates_segment()	n_jobs=1	Runs k-means clustering iterations in parallel via joblib.Parallel

Internal optimization

• _hrv_dfa() in hrv_nonlinear.py: monofractal and multifractal DFA now run concurrently via concurrent.futures.ThreadPoolExecutor (always-on, with sequential fallback on failure).

Other additions

• tests/tests_parallel.py: 4 new tests verifying parallel results match sequential results for bio_process, ecg_findpeaks ProMAC, and entropy_multiscale.
• benchmarks/bench_parallel.py: Benchmark script comparing sequential vs parallel performance.
• README.rst: Added "Parallel Processing" section with usage examples.

Benchmark Results

Measured on a 12-core machine (AMD/Intel), Python 3.12, averaged over 3 runs:

============================================================
  SUMMARY
============================================================
  Benchmark                                      Seq (s)  Par (s)  Speedup
  --------------------------------------------- -------- -------- --------
  bio_process (5 modalities, 5min @ 1000Hz)        7.669    4.994    1.54x
  ecg_findpeaks ProMAC (5min @ 1000Hz)             6.582    4.579    1.44x
  entropy_multiscale (30s @ 500Hz, 30 scales)      0.281    0.133    2.11x
  HRV nonlinear (DFA threads)                      0.636      -        -
============================================================

Speedups scale with workload size; longer recordings and more channels yield larger gains.

Design Decisions

• n_jobs=1 as default everywhere to preserve existing behavior and avoid surprises.
• bio_process uses parallel=False (bool) instead of n_jobs because it spawns one process per signal type (up to 6), not an arbitrary worker pool.
• joblib is used for the n_jobs pattern (consistent with scikit-learn convention). concurrent.futures is used where stdlib is sufficient (bio_process, HRV DFA).
• No new hard dependencies added.

Checklist

• I have read the CONTRIBUTING file.
• My PR is targeted at the dev branch (and not towards the master branch).
• I ran the CODE CHECKS on the files I added or modified and fixed the errors.
• I have added the newly added features to News.rst (if applicable)

[gemini-code-assist]

Code Review

This pull request introduces optional parallel processing support across several computationally intensive functions in NeuroKit2, including bio_process, ecg_findpeaks (ProMAC), entropy_multiscale, eeg_power, and microstates_segment, as well as internal thread-based parallelism for DFA in hrv_nonlinear. The review feedback highlights several areas for improvement: ensuring robust error handling in the parallel ProMAC implementation to prevent total failure on individual method errors, optimizing resource allocation in ProcessPoolExecutor by limiting workers to the task count, reducing overhead by batching pandas concatenation operations, and refactoring duplicated logic between sequential and parallel execution paths to improve maintainability and performance.

[gemini-code-assist]

The parallel implementation of ProMAC does not handle exceptions for individual methods, unlike the sequential version. If one peak detection method fails, the entire joblib.Parallel call will raise an exception, causing ecg_findpeaks to fail. In the sequential version, errors are caught and stored in error_list, allowing the process to continue with other methods. You should wrap the logic inside the worker in a try...except block to maintain this robustness.

[gemini-code-assist]

When using ProcessPoolExecutor without specifying max_workers, it defaults to the number of processors on the machine. Since the number of tasks is fixed and small (at most 6), spawning a large pool is inefficient. It is recommended to limit the number of workers to the number of tasks to reduce overhead.

Suggested change

	with concurrent.futures.ProcessPoolExecutor() as executor:
	with concurrent.futures.ProcessPoolExecutor(max_workers=len(_tasks)) as executor:

[gemini-code-assist]

Repeatedly calling pd.concat in a loop is computationally expensive because it involves copying the entire DataFrame at each step. A more efficient approach is to collect all DataFrames in a list and perform a single concatenation after the loop.

    # Assemble results in consistent order
    to_concat = [bio_df]
    for name, _, _ in _tasks:
        signals_df, info_dict = results[name]
        bio_info.update(info_dict)
        to_concat.append(signals_df)
    bio_df = pd.concat(to_concat, axis=1)

[gemini-code-assist]

The sequential and parallel execution paths contain duplicated logic for calling _entropy_multiscale. This redundancy increases the maintenance burden. You can refactor this by defining a helper function or a task generator and then choosing the execution method based on n_jobs.

    def _run(scale_factor):
        return _entropy_multiscale(
            signal,
            scale=scale_factor,
            coarsegraining=coarsegraining,
            algorithm=algorithm,
            dimension=dimension,
            tolerance=info["Tolerance"],
            refined=refined,
            **kwargs,
        )

    if n_jobs == 1:
        # Sequential execution (original behavior)
        info["Value"] = np.array([_run(s) for s in info["Scale"]])
    else:
        # Parallel execution via joblib
        try:
            import joblib
        except ImportError as e:
            raise ImportError(
                "NeuroKit error: entropy_multiscale(): the 'joblib' module is required "
                "for parallel execution. Please install it first (`pip install joblib`).",
            ) from e

        info["Value"] = np.array(
            joblib.Parallel(n_jobs=n_jobs)(
                joblib.delayed(_run)(s) for s in info["Scale"]
            )
        )

[gemini-code-assist]

The _run_method helper duplicates the convolution and Gaussian kernel logic found in _ecg_findpeaks_promac_addconvolve. Additionally, the Gaussian shape is recalculated for every method, even though it only depends on the sampling rate and gaussian_sd. It is better to calculate the kernel once and reuse the existing convolution logic to ensure consistency.

[feruzoripov]

all comments were addressed in new commits

[codecov-commenter]

Codecov Report

❌ Patch coverage is 83.15789% with 16 lines in your changes missing coverage. Please review.
✅ Project coverage is 59.02%. Comparing base (608773f) to head (7c1ca9b).
⚠️ Report is 6 commits behind head on dev.

Files with missing lines	Patch %	Lines
neurokit2/microstates/microstates_segment.py	56.52%	10 Missing ⚠️
neurokit2/bio/bio_process.py	90.62%	3 Missing ⚠️
neurokit2/ecg/ecg_findpeaks.py	87.50%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##              dev    #1187      +/-   ##
==========================================
+ Coverage   57.85%   59.02%   +1.16%     
==========================================
  Files         310      310              
  Lines       15701    15757      +56     
==========================================
+ Hits         9084     9300     +216     
+ Misses       6617     6457     -160     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[feruzoripov]

Should I request a review from someone?

[DerAndereJohannes]

Hi feruzoripov, thank you for the detailed PR with benchmarks for your changes! Also excuse the late response. This is very beneficial for the embarrassingly parallel sections of the code like in bio_process. I would like to preface my review with saying that I do think that the default should always be sequential regardless of any speedups mostly due to keeping the default as stable as possible and to prevent overhead if people put in small signals (Which would make it slower!).

Now a few things that are open to discussion from my side:

1. parallel_run.py: NeuroKit already has a parallel_run() utility in neurokit2/misc/ that wraps joblib.Parallel. For the n_jobs additions (ecg_findpeaks, entropy_multiscale, eeg_power, microstates_segment), was there a reason not to use that directly rather than reimplementing the same joblib.Parallel boilerplate in each function? Where it fits naturally (e.g., looping over scale factors in entropy_multiscale, channels in eeg_power, or methods in ProMAC), I would personally prefer to see parallel_run() used directly.

2. n_jobs in ecg_findpeaks: Adding n_jobs as a top-level parameter here is misleading since it only
   does anything when method="promac". For every other method it would be silently ignored. I think a cleaner approach would be to let it pass through **kwargs down to _ecg_findpeaks_promac(), which can consume it with n_jobs = kwargs.pop("n_jobs", 1) before forwarding the remaining kwargs to the individual peak detection methods. That way the public signature stays clean and there's no risk of users passing n_jobs with a non-ProMAC method and expecting something to happen.

3. bio_process and _hrv_dfa backends: The PR introduces concurrent.futures.ProcessPoolExecutor for bio_process and ThreadPoolExecutor for _hrv_dfa, alongside joblib everywhere else. Since joblib already supports both workload types (backend="loky" for processes, backend="threading" for threads), I would much prefer to standardize on joblib throughout. The practical reason for bio_process specifically: I don't think ProcessPoolExecutor detects nested parallelism contexts, so if a user runs parallel_run over subjects and each subject's bio_process also parallelizes internally, you risk spawning processes inside workers with no coordination.

4. _hrv_dfa threading: The broad silent except Exception fallback is always a red flag for me. Any real error inside fractal_dfa would be swallowed without the user knowing. Beyond that, only 2 tasks are being threaded, so the overhead of creating a ThreadPoolExecutor (or joblib as per the comment 3) could easily outweigh the benefit, especially for shorter recordings. Whether threading helps at all depends on whether fractal_dfa releases the GIL, which isn't obvious. Could you add a proper before/after benchmark (I see most of the results here are -)? Similarly, I would suggest making this opt-in via n_jobs (if the benchmark shows that it is worth it) and dropping the silent fallback.

Again, thank you for taking the time for your contribution! I look forward to your responses.

I would also be happy if others posted their opinions.

[feruzoripov]

Thanks for the thorough review @DerAndereJohannes! All four points addressed:

1. Use parallel_run(): Agreed, replaced all raw joblib.Parallel boilerplate and concurrent.futures usage with parallel_run() across bio_process, ecg_findpeaks ProMAC, entropy_multiscale, eeg_power, and microstates_segment. Single utility, consistent behavior.

2. n_jobs in ecg_findpeaks: Removed from the public signature. It now flows through **kwargs and is consumed inside _ecg_findpeaks_promac() via kwargs.pop("n_jobs", 1), so it's invisible to non-ProMAC methods.

3. Standardize on joblib: Dropped oncurrent.futures entirely. bio_process now uses parallel_run() (joblib backend), which also means nested parallelism is handled by joblib's built-in coordination (loky).

4. _hrv_dfa threading: Reverted to the original sequential code. You're right, the silent except Exception was a red flag, and the benefit of threading 2 tasks was questionable without a clear GIL-release guarantee. Removed entirely rather than making it opt-in, since the benchmark didn't show meaningful gains.

And yes, fully agree that defaults should always be sequential, that's been the case throughout (n_jobs=1 / parallel=False).