Add possibility to fill missing parcels when extracting timeseries with "fill_missing" flag/argument and bump supported python version to 3.10-3.14

.circleci/config.yml

@@ -10,19 +10,19 @@ orbs:
   # so you dont have to copy and paste it everywhere.
   # See the orb documentation here: https://circleci.com/developer/orbs/orb/circleci/python
   python: circleci/python@2.0.3
-  codecov: codecov/codecov@3.2.2
+  codecov: codecov/codecov@6.1.0
 
 # Define a job to be invoked later in a workflow.
 # See: https://circleci.com/docs/2.0/configuration-reference/#jobs
 jobs:
-  test39: # This is the name of the job, feel free to change it to better match what you're trying to do!
+  test311: # This is the name of the job, feel free to change it to better match what you're trying to do!
     # These next lines defines a Docker executors: https://circleci.com/docs/2.0/executor-types/
     # You can specify an image from Dockerhub or use one of the convenience images from CircleCI's Developer Hub
     # A list of available CircleCI Docker convenience images are available here: https://circleci.com/developer/images/image/cimg/python
     # The executor is the environment in which the steps below will be executed - below will use a python 3.6.14 container
     # Change the version below to your required version of python
     docker:
-      - image: cimg/python:3.9
+      - image: cimg/python:3.11
     working_directory: /tmp/src/nigsp
     resource_class: medium
     # Checkout the code as the first step. This is a dedicated CircleCI step.
@@ -46,7 +46,7 @@ jobs:
           command: |
             pytest --cov=./nigsp
             mkdir /tmp/src/coverage
-            mv ./.coverage /tmp/src/coverage/.coverage.py37
+            mv ./.coverage /tmp/src/coverage/.coverage.py311
       - store_artifacts:
           path: /tmp/src/coverage
       # Persist the specified paths (workspace/echo-output) into the workspace for use in downstream job.
@@ -56,11 +56,11 @@ jobs:
           root: /tmp
           # Must be relative path from root
           paths:
-            - src/coverage/.coverage.py37
+            - src/coverage/.coverage.py311
 
-  test312:
+  test314:
     docker:
-      - image: cimg/python:3.12
+      - image: cimg/python:3.14
     working_directory: /tmp/src/nigsp
     resource_class: medium
     steps:
@@ -75,13 +75,13 @@ jobs:
           command: |
             pytest --cov=./nigsp
             mkdir /tmp/src/coverage
-            mv ./.coverage /tmp/src/coverage/.coverage.py310
+            mv ./.coverage /tmp/src/coverage/.coverage.py314
       - store_artifacts:
           path: /tmp/src/coverage
       - persist_to_workspace:
           root: /tmp
           paths:
-            - src/coverage/.coverage.py310
+            - src/coverage/.coverage.py314
 
   merge_coverage:
     working_directory: /tmp/src/nigsp
@@ -95,27 +95,29 @@ jobs:
       - python/install-packages:
           args: coverage
           pkg-manager: pip
+      - run:
+          name: Install gpg (needed for codecov)
+          command: sudo apt-get update && sudo apt-get install -y gpg curl
       - run:
           name: Merge coverage files
           command: |
-            sudo apt update && sudo apt install curl
             cd /tmp/src/coverage/
             coverage combine
             coverage xml
       - store_artifacts:
           path: /tmp/src/coverage
       - codecov/upload:
-          file: /tmp/src/coverage/coverage.xml
+          files: /tmp/src/coverage/coverage.xml
 
 # Invoke jobs via workflows
 # See: https://circleci.com/docs/2.0/configuration-reference/#workflows
 workflows:
   build_test: # This is the name of the workflow, feel free to change it to better match your workflow.
     # Inside the workflow, you define the jobs you want to run.
     jobs:
-      - test39
-      - test312
+      - test311
+      - test314
       - merge_coverage:
           requires:
-            - test39
-            - test312
+            - test311
+            - test314

nigsp/blocks.py

@@ -16,7 +16,7 @@
 LGR = logging.getLogger(__name__)
 
 
-def nifti_to_timeseries(fname, atlasname):
+def nifti_to_timeseries(fname, atlasname, fill_missing=False):
     """Read a nifti file and returns a normalised timeseries from an atlas.
 
     Parameters
@@ -25,6 +25,8 @@ def nifti_to_timeseries(fname, atlasname):
         Filename (and path) of a functional timeseries nifti dataset.
     atlasname : str | os.PathLike
         Filename (and path) of an atlas nifti dataset.
+    fill_missing : bool, optional
+        If True, fill missing atlas parcels with timeseries of 0. Default is False.
 
     Returns
     -------
@@ -42,7 +44,7 @@ def nifti_to_timeseries(fname, atlasname):
     atlas, amask, img = io.load_nifti_get_mask(atlasname, is_mask=True, ndim=3)
     mask *= amask
 
-    timeseries = nifti.apply_atlas(data, atlas, mask)
+    timeseries = nifti.apply_atlas(data, atlas, mask, fill_missing)
 
     return timeseries, atlas, img
 

nigsp/cli/run.py

@@ -119,6 +119,15 @@ def _get_parser():
         ),
         default="median",
     )
+    opt_proc.add_argument(
+        "-fill--fill_missing",
+        dest="fill_missing",
+        action="store_true",
+        help=(
+            "Fill missing atlas parcels with timeseries of 0s when extracting "
+            "timeseries. This will not check for a last parcel missing."
+        ),
+    )
 
     opt_out = parser.add_argument_group("Optional Arguments for output")
     opt_out.add_argument(

nigsp/conftest.py

@@ -1,5 +1,6 @@
 from __future__ import annotations  # c.f. PEP 563, PEP 649
 
+import logging
 import os
 import ssl
 from typing import TYPE_CHECKING
@@ -53,10 +54,8 @@ def fetch_file(osf_id, path, filename):
     # this three lines make tests downloads work in windows
     if os.name == "nt":
         orig_sslsocket_init = ssl.SSLSocket.__init__
-        ssl.SSLSocket.__init__ = (
-            lambda *args, cert_reqs=ssl.CERT_NONE, **kwargs: orig_sslsocket_init(
-                *args, cert_reqs=ssl.CERT_NONE, **kwargs
-            )
+        ssl.SSLSocket.__init__ = lambda *args, cert_reqs=ssl.CERT_NONE, **kwargs: (
+            orig_sslsocket_init(*args, cert_reqs=ssl.CERT_NONE, **kwargs)
         )
         ssl._create_default_https_context = ssl._create_unverified_context
     url = f"https://osf.io/{osf_id}/download"
@@ -104,3 +103,19 @@ def sc_mtx(testdir):
 @fixture(scope="function")
 def timeseries(testdir):
     return fetch_file("ay8df", testdir, "func.mat")
+
+
+@fixture(autouse=True)
+def cleanup_logging_handlers():
+    """Close and remove all logging handlers after each test run."""
+    yield
+
+    # Target both the root logger and any package-specific loggers
+    for name in [None, "nigsp"]:
+        logger = logging.getLogger(name)
+        for handler in logger.handlers[:]:
+            handler.close()
+            logger.removeHandler(handler)
+
+    # Flush and close all remaining global logging streams
+    logging.shutdown()

nigsp/operations/laplacian.py

@@ -58,8 +58,8 @@ def compute_laplacian(mtx, negval="absolute", selfloops=False):
     -----
     https://en.wikipedia.org/wiki/Laplacian_matrix
     """
-    mtx = deepcopy(mtx)
     if mtx.min() < 0:
+        mtx = deepcopy(mtx)
         if negval == "absolute":
             mtx = abs(mtx)
         elif negval == "remove":

nigsp/operations/nifti.py

@@ -165,7 +165,7 @@ def unmask(data, mask, shape=None, asdata=None):
     return out
 
 
-def apply_atlas(data, atlas, mask=None):
+def apply_atlas(data, atlas, mask=None, fill_missing=False):
     """
     Extract average timeseries from an atlas.
 
@@ -179,6 +179,8 @@ def apply_atlas(data, atlas, mask=None):
     mask : None or numpy.ndarray, optional
         A 2- or 3- D matrix representing a mask, all voxels == 0 are excluded from the
         computation.
+    fill_missing : bool, optional
+        If True, fill missing atlas parcels with timeseries of 0. Default is False.
 
     Returns
     -------
@@ -228,11 +230,24 @@ def apply_atlas(data, atlas, mask=None):
     labels = labels[labels > 0]
     LGR.info(f"Labels: {labels}, numbers: {len(labels)}")
     # Initialise dataframe and dictionary for series
-    parcels = np.empty([len(labels), data.shape[-1]], dtype="float32")
+
+    if fill_missing:
+        all_labels = np.arange(1, labels[-1] + 1)
+
+        if all_labels.size != labels.size:
+            LGR.warning(
+                f"Found {labels.size} labels, expected {all_labels.size}. Filling "
+                "missing timeseries."
+            )
+            labels = all_labels
+
+    parcels = np.zeros([len(labels), data.shape[-1]], dtype="float32")
 
     # Compute averages
     for n, label in enumerate(labels):
-        parcels[n, :] = data[atlas == label].mean(axis=0)
+        ts_mtx = data[atlas == label]
+        if ts_mtx.size > 0:
+            parcels[n, :] = ts_mtx.mean(axis=0)
 
     return parcels
 

nigsp/operations/surrogates.py

@@ -65,12 +65,12 @@ def random_sign(eigenvec, n_surr=1000, seed=42, stack=False):
         if eigenvec.ndim < 3:
             r_sign = rng.integers(0, 1, eigenvec.shape[0], endpoint=True)
             r_sign[r_sign == 0] = -1
-            rand_evec[..., i] = eigenvec * r_sign
+            rand_evec[..., i] = (eigenvec * r_sign).real
         else:
             for j in range(eigenvec.shape[2]):
                 r_sign = rng.integers(0, 1, eigenvec.shape[0], endpoint=True)
                 r_sign[r_sign == 0] = -1
-                rand_evec[:, :, j, i] = eigenvec[..., j] * r_sign
+                rand_evec[:, :, j, i] = (eigenvec[..., j] * r_sign).real
 
     if stack:
         rand_evec = np.append(rand_evec, eigenvec[..., np.newaxis], axis=-1)

nigsp/operations/timeseries.py

@@ -264,11 +264,11 @@ def graph_fourier_transform(timeseries, eigenvec, energy=False, mean=False):
     """
     timeseries = timeseries.squeeze()
     if timeseries.ndim < 3:
-        proj = eigenvec.conj().T @ timeseries
+        proj = (eigenvec.conj().T @ timeseries).real
     else:
         temp_ts, proj = prepare_ndim_iteration(timeseries, 2)
         for i in range(temp_ts.shape[-1]):
-            proj[:, :, i] = eigenvec.conj().T @ np.squeeze(temp_ts[:, :, i])
+            proj[:, :, i] = (eigenvec.conj().T @ np.squeeze(temp_ts[:, :, i])).real
         if timeseries.ndim > 3:
             proj = proj.reshape(timeseries.shape)
 

nigsp/tests/test_nifti.py

@@ -99,6 +99,24 @@ def test_apply_atlas():
     assert rm.shape == cm.shape
     assert (rm == cm).all()
 
+    a = asarray([1, 2, 2, 4, 4, 0])
+
+    r = nifti.apply_atlas(d, a, fill_missing=False)
+
+    assert r.shape == c.shape
+    assert (r == c).all()
+
+    c = zeros((4, 10), dtype="float32")
+
+    c[0, :] = d[0, :]
+    c[1, :] = d[1:3, :].mean(axis=0)
+    c[3, :] = d[3:5, :].mean(axis=0)
+
+    r = nifti.apply_atlas(d, a, fill_missing=True)
+
+    assert r.shape == c.shape
+    assert (r == c).all()
+
 
 def test_unfold_atlas():
     m = asarray([0, 1, 1, 1, 1, 1])

nigsp/tests/test_surrogates.py

@@ -22,7 +22,7 @@ def test_random_sign():
     for i in range(n_surr):
         r_sign = rng.integers(0, 1, eigenvec.shape[0], endpoint=True)
         r_sign[r_sign == 0] = -1
-        rand_evec[..., i] = eigenvec * r_sign
+        rand_evec[..., i] = (eigenvec * r_sign).real
     rand_evec = np.append(rand_evec, eigenvec[..., np.newaxis], axis=-1)
 
     assert rs.ndim == 2

nigsp/workflow.py

@@ -75,6 +75,7 @@ def nigsp(
     outname=None,
     outdir=None,
     comp_metric=[],  # noqa: B006
+    fill_missing=False,
     index="median",
     surr_type=None,
     n_surr=1000,
@@ -114,6 +115,8 @@ def nigsp(
     comp_metric : list or None, optional
         List of metrics that should be computed. If empty (default), compute all
         metrics available.
+    fill_missing : bool, optional
+        If True, fill missing atlas parcels with timeseries of 0. Default is False.
     index : 'median' or int, optional
         The index of the eigenvector/harmonic of the graph to split the graph in
         multiple parts, or the method to find this index. Currently supports
@@ -334,7 +337,9 @@ def nigsp(
     timeseries = []
     for f in fname:
         if func_is["nifti"] and atlas_is["nifti"]:
-            t, atlas, img = blocks.nifti_to_timeseries(f, atlasname)
+            t, atlas, img = blocks.nifti_to_timeseries(
+                f, atlasname, fill_missing=fill_missing
+            )
         elif func_is["nifti"] and atlas_is["nifti"] is False:
             raise NotImplementedError(
                 "To work with functional file(s) of nifti format, "
@@ -470,8 +475,10 @@ def nigsp(
                 pass
 
     LGR.info(f"End of workflow, find results in {outdir}.")
-    LGR.removeHandler(log_handler)
     log_handler.close()
+    LGR.removeHandler(log_handler)
+    sh.close()
+    LGR.removeHandler(sh)
 
     return 0
 

setup.cfg

@@ -9,8 +9,6 @@ classifiers =
     Development Status :: 4 - Beta
     Intended Audience :: Science/Research
     License :: OSI Approved :: Apache Software License
-    Programming Language :: Python :: 3.8
-    Programming Language :: Python :: 3.9
     Programming Language :: Python :: 3.10
     Programming Language :: Python :: 3.11
     Programming Language :: Python :: 3.12
@@ -42,7 +40,8 @@ mat =
 nifti =
     nibabel
 viz =
-    matplotlib>=3.1.1
+    matplotlib>=3.10.0
+    pyparsing>=3.0.0
     nilearn>=0.7.0
 all =
     %(nifti)s