Merge pull request #77 from MICA-MNI/vtk9.2

Vtk9.2

.github/workflows/tagged_release.yml

@@ -0,0 +1,45 @@
+name: "tagged_release"
+
+on:
+  push:
+    branches:
+      - master
+    tags:
+      - "v*"
+
+jobs:
+  tagged-release:
+    name: "Tagged Release"
+    runs-on: "ubuntu-latest"
+
+    steps:
+    - uses: actions/checkout@v2
+
+    - name: Set up Python 3.8.
+      uses: actions/setup-python@v2
+      with:
+        python-version: 3.8
+
+    - name: Install Python BrainSpace.
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install --upgrade setuptools wheel
+
+    - name: Build binary wheel and tarball.
+      run: |
+        python setup.py sdist bdist_wheel
+
+    - name: Create Github release.
+      uses: "marvinpinto/action-automatic-releases@latest"
+      with:
+        repo_token: "${{ secrets.GITHUB_TOKEN }}"
+        prerelease: false
+        files: |
+          LICENSE
+          dist/*
+
+    - name: Publish to PyPi.
+      uses: pypa/gh-action-pypi-publish@release/v1
+      with:
+        password: ${{ secrets.PYPI_TOKEN }}
+        verbose: true

CHANGELOG.txt

@@ -1,11 +1,16 @@
-v0.1.4, April 2022 -- Bug Fixes
+v0.1.4, August 2022 -- Support for vtk 9.2
 MATLAB:
   - Fixed the uploading of tutorial data to FileExchange.
-  - Fixed a bug that caused the detection of connected affinity matrices to error.
+  - Fixed a bug that caused the detection of connected affinity matrices to
+  error.
 
 Python:
-  - GradientMaps.fit() now returns an error when the affinity matrix contains NaN or Inf.
-  - Added a warning when plot_hemispheres() is run without a display. 
+  - GradientMaps.fit() now returns an error when the affinity matrix contains
+  NaN or Inf.
+  - Added a warning when plot_hemispheres() is run without a display.
+  - Added support for vtk 9.2
+  - For plotting, array_name=None now plots surfaces without any array data.
+
 
 v0.1.3, January 2022 -- Fix dependency error
 Python:

brainspace/_version.py

@@ -1,3 +1,3 @@
 """BrainSpace version"""
 
-__version__ = '0.1.3'
+__version__ = '0.1.4'

brainspace/gradient/gradient.py

@@ -29,7 +29,7 @@ def _fit_one(x, app, kernel, n_components, random_state, gamma=None,
         Inverse kernel width. Only used if ``kernel`` == 'gaussian'.
         If None, ``gamma=1/n_feat``. Default is None.
     sparsity : float, optional
-        Proportion of smallest elements to zero-out for each row.
+        Proportion of the smallest elements to zero-out for each row.
         Default is 0.9.
     kwargs : kwds, optional
         Additional keyword parameters passed to the embedding approach.
@@ -45,9 +45,9 @@ def _fit_one(x, app, kernel, n_components, random_state, gamma=None,
     a = compute_affinity(x, kernel=kernel, sparsity=sparsity, gamma=gamma)
 
     if np.isnan(a).any() or np.isinf(a).any():
-        raise ValueError(
-            "Affinity matrix contains NaN or Inf values. Common causes of this include NaNs/Infs or rows of zeros in the input matrix."
-        )
+        raise ValueError('Affinity matrix contains NaN or Inf values. Common '
+                         'causes of this include NaNs/Infs or rows of zeros '
+                         'in the input matrix.')
 
     kwds_emb = {'n_components': n_components, 'random_state': random_state}
     kwds_emb.update(kwargs)
@@ -87,7 +87,7 @@ class GradientMaps(BaseEstimator):
         If None, use input matrix. Default is None.
     alignment : {'procrustes', 'joint'}, object or None
         Alignment approach. Only used when two or more datasets are provided.
-        If None, no alignment is peformed. If `object`, it accepts an instance
+        If None, no alignment is performed. If `object`, it accepts an instance
         of :class:`.ProcrustesAlignment`. Default is None.
 
         - If 'procrustes', datasets are aligned using generalized procrustes
@@ -134,7 +134,7 @@ def fit(self, x, gamma=None, sparsity=0.9, n_iter=10, reference=None,
             Inverse kernel width. Only used if ``kernel == 'gaussian'``.
             If None, ``gamma=1/n_feat``. Default is None.
         sparsity : float, optional
-            Proportion of smallest elements to zero-out for each row.
+            Proportion of the smallest elements to zero-out for each row.
             Default is 0.9.
         n_iter : int, optional
             Number of iterations for procrustes alignment. Default is 10.

brainspace/gradient/utils.py

@@ -38,8 +38,8 @@ def is_symmetric(x, tol=1E-10):
     if ssp.issparse(x):
         if x.format not in ['csr', 'csc', 'coo']:
             x = x.tocoo(copy=False)
-        dif1 = x - x.T
-        return np.all(np.abs(dif1.data) < tol)
+        dif = x - x.T
+        return np.all(np.abs(dif.data) < tol)
 
     return np.allclose(x, x.T, atol=tol)
 
@@ -141,7 +141,7 @@ def _dominant_set_dense(s, k, is_thresh=False, norm=False, copy=True):
 
 
 def dominant_set(s, k, is_thresh=False, norm=False, copy=True, as_sparse=True):
-    """Keep largest elements for each row. Zero-out the rest.
+    """Keep the largest elements for each row. Zero-out the rest.
 
     Parameters
     ----------
@@ -183,3 +183,94 @@ def dominant_set(s, k, is_thresh=False, norm=False, copy=True, as_sparse=True):
         return _dominant_set_sparse(s, k, is_thresh=is_thresh, norm=norm)
 
     return _dominant_set_dense(s, k, is_thresh=is_thresh, norm=norm, copy=copy)
+
+
+def ravel_symmetric(x, with_diagonal=False):
+    """Return the flattened upper triangular part of a symmetric matrix.
+
+    Parameters
+    ----------
+    x : 2D ndarray or sparse matrix, shape=(n, n)
+        Input array.
+    with_diagonal : bool, optional
+        If True, also return diagonal elements. Default is False.
+
+    Returns
+    -------
+    output : 1D ndarray, shape (n_feat,)
+        The flattened upper triangular part of `x`. If with_diagonal
+        is True, ``n_feat = n * (n + 1) / 2`` and
+        ``n_feat = n * (n - 1) / 2`` otherwise.
+
+    """
+
+    n = x.shape[0]
+    k = 0 if with_diagonal else -1
+    mask_lt = np.tri(n, k=k, dtype=np.bool)
+
+    if ssp.issparse(x) and not ssp.isspmatrix_csc(x):
+        x = x.tocsc(copy=False)
+
+    return x[mask_lt.T]
+
+
+def unravel_symmetric(x, size, as_sparse=False, part='both', fmt='csr'):
+    """Build symmetric matrix from array with upper triangular elements.
+
+    Parameters
+    ----------
+    x : 1D ndarray
+        Input data with elements to go in the upper triangular part.
+    size : int
+        Number of rows/columns of matrix.
+    as_sparse : bool, optional
+        Return a sparse matrix. Default is False.
+    part: {'both', 'upper', 'lower'}, optional
+        Build matrix with elements if both or just on triangular part.
+        Default is both.
+    fmt: str, optional
+        Format of sparse matrix. Only used if ``as_sparse=True``.
+        Default is 'csr'.
+
+    Returns
+    -------
+    sym : 2D ndarray or sparse matrix, shape = (size, size)
+        Array with the lower/upper or both (symmetric) triangular parts
+        built from `x`.
+
+    """
+
+    k = 1
+    if (size * (size + 1) // 2) == x.size:
+        k = 0
+    elif (size * (size - 1) // 2) != x.size:
+        raise ValueError('Cannot unravel data. Wrong size.')
+
+    shape = (size, size)
+    if as_sparse:
+        mask = x != 0
+        x = x[mask]
+
+        idx = np.triu_indices(size, k=k)
+        idx = [idx1[mask] for idx1 in idx]
+        if part == 'lower':
+            idx = idx[::-1]
+        elif part == 'both':
+            idx = np.concatenate(idx), np.concatenate(idx[::-1])
+            x = np.tile(x, 2)
+
+        xs = ssp.coo_matrix((x, idx), shape=shape)
+        if fmt != 'coo':
+            xs = xs.asformat(fmt, copy=False)
+
+    else:
+        mask_lt = np.tri(size, k=-k, dtype=np.bool)
+        xs = np.zeros(shape, dtype=x.dtype)
+
+        xs[mask_lt.T] = x
+        if part == 'both':
+            xs.T[mask_lt.T] = x
+        elif part == 'lower':
+            xs = xs.T
+
+    return xs

brainspace/mesh/array_operations.py

@@ -10,7 +10,7 @@
 
 import numpy as np
 from scipy.stats import mode
-from scipy.spatial import cKDTree
+from scipy.spatial import KDTree
 from scipy.sparse.csgraph import laplacian, connected_components
 
 from sklearn.utils.extmath import weighted_mode
@@ -624,31 +624,41 @@ def propagate_labeling(surf, labeling, no_label=np.nan, mask=None, alpha=0.99,
     # Graph matrix
     if mode == 'connectivity':
         adj = me.get_ring_adjacency(surf, n_ring=n_ring, include_self=False,
-                                    dtype=np.float)
+                                    dtype=np.float, mask=mask)
     else:
-        adj = me.get_ring_distance(surf, n_ring=n_ring, dtype=np.float)
+        adj = me.get_ring_distance(surf, n_ring=n_ring, dtype=np.float,
+                                   mask=mask)
         adj.data[:] = np.exp(-adj.data/n_ring**2)
 
     if mask is not None:
         adj = adj[mask][:, mask]
 
-    graph_matrix = -alpha * laplacian(adj, normed=True)
-    diag_mask = (graph_matrix.row == graph_matrix.col)
-    graph_matrix.data[diag_mask] = 0.0
+    # graph_matrix = -alpha * laplacian(adj, normed=True)
+    graph = laplacian(adj, normed=True)
+    graph.data *= -alpha
+    graph.data[graph.row == graph.col] = 0
+    graph.eliminate_zeros()
+    graph = graph.tocsr(copy=False)
+    # diag_mask = (graph_matrix.row == graph_matrix.col)
+    # graph_matrix.data[diag_mask] = 0.0
 
     # Label distributions and label static
     lab_dist = np.zeros((n_pts, n_labs))
     lab_dist[np.argwhere(labeled)[:, 0], idx_lab] = 1
 
-    lab_static = lab_dist.copy()
-    lab_static *= 1 - alpha
+    # lab_static = lab_dist.copy()
+    # lab_static *= 1 - alpha
+    lab_static = (1 - alpha) * lab_dist
 
     # propagation
     lab_dist_perv = lab_dist
     for i in range(n_iter):
-        lab_dist = graph_matrix.dot(lab_dist) + lab_static
+        lab_dist = graph.dot(lab_dist)
+        lab_dist += lab_static
+        # lab_dist = graph_matrix.dot(lab_dist) + lab_static
 
-        if np.linalg.norm(lab_dist - lab_dist_perv, 'fro') < tol:
+        lab_dist_perv -= lab_dist
+        if np.linalg.norm(lab_dist_perv, 'fro') < tol:
             break
 
         lab_dist_perv = lab_dist
@@ -823,15 +833,14 @@ def _get_pids_sphere(source, target, source_mask=None, target_mask=None):
 
 
 def _get_pids_naive(source, target, k=1, source_mask=None, target_mask=None,
-                    return_weights=True, n_jobs=1):
-    """Resampling based on k nearest points."""
+                    return_weights=True):
+    """Resampling based on the k nearest points."""
 
     sp = me.get_points(source, mask=source_mask)
     tp = me.get_points(target, mask=target_mask)
 
-    tree = cKDTree(sp, leafsize=20, compact_nodes=False, copy_data=False,
-                   balanced_tree=False)
-    dist, pids = tree.query(tp, k=k, eps=0, n_jobs=n_jobs)
+    tree = KDTree(sp, leafsize=20)
+    dist, pids = tree.query(tp, k=k, eps=0)
 
     if return_weights:
         return pids, 1 / dist
@@ -899,6 +908,9 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
     """
     opt = ['mean', 'mode', 'weighted_mean', 'weighted_mode']
 
+    if n_jobs != 1:
+        warnings.warn('The n_jobs parameter is deprecated and will be removed '
+                      'in a future version', DeprecationWarning)
     is_list = True
     if not isinstance(data, list):
         data = [data]
@@ -918,7 +930,7 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
             use_weights = True
 
         pids = _get_pids_naive(source, target, k=k, source_mask=source_mask,
-                               target_mask=target_mask, n_jobs=n_jobs,
+                               target_mask=target_mask,
                                return_weights=use_weights)
         if use_weights:
             pids, w = pids
@@ -942,7 +954,7 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
         if k == 1:
             feat = d[pids]
         elif red_func[i] == 'mean':
-            feat = np.mean(d[pids], axis=1)
+            feat = np.nanmean(d[pids], axis=1)
         elif red_func[i] == 'weighted_mean':
             feat = np.average(d[pids], weights=w, axis=1)
         elif red_func[i] == 'mode':
@@ -954,7 +966,7 @@ def resample_pointdata(source, target, data, is_sphere=False, source_mask=None,
             raise ValueError('Unknown red_func: {0}'.format(red_func[i]))
 
         if target_mask is not None:
-            feat = map_to_mask(feat, mask=target_mask, fill=fill)
+            feat = map_to_mask(feat, mask=target_mask, axis=1, fill=fill)
         resampled[i] = feat
 
     if append and key is not None:

brainspace/mesh/mesh_cluster.py

@@ -90,7 +90,8 @@ def cluster_points(surf, n_clusters=100, is_size=False, mask=None,
     # Find clusters
     if approach == 'kmeans':
         if n_jobs != 1:
-            warn("The n_jobs parameter is deprecated and will be removed in a future version", DeprecationWarning)
+            warn('The n_jobs parameter is deprecated and will be removed '
+                 'in a future version', DeprecationWarning)
         _, cluster_labs, _ = k_means(evs, n_clusters=n_clusters,
                                      random_state=random_state, n_init=n_init)
     else: