Merge pull request #9 from SNG-Newy/8-pypi-install

pypi upload

MANIFEST.in

@@ -1,3 +1,4 @@
 include README.rst LICENSE
 recursive-include eigenstrapping/datasets *
 include eigenstrapping/datasets/osf.json
+exclude *egg-info* dist __pycache__ 

README.rst

@@ -18,7 +18,7 @@
 The ``eigenstrapping`` toolbox is designed to help researchers generate
 statistically-rigorous models for null hypothesis testing between brain maps using
 non-local spectral shape descriptors - or geometric eigenmodes.
-Documentation can be found `here <https://eigenstrapping.readthedocs.io/en/latest/>`_. Read our preprint `here <https://dx.doi.org/10.1101/2024.02.07.579070>`_!
+Documentation can be found `here <https://eigenstrapping.readthedocs.io/en/latest/>`_. Read our preprint `<https://dx.doi.org/10.1101/2024.02.07.579070>`_
 
 .. image:: ./docs/_static/main_figure.jpg
   :scale: 20%
@@ -62,7 +62,7 @@ To run eigenstrapping, the following Python packages are required (these should
 * `netneurotools <https://github.com/netneurolab/netneurotools>`_
 
 ``nibabel`` and ``nilearn`` are required for surfaces and volumes. ``matplotlib``
-is only required for fitting plots in :mod:`eigenstrapping.fit` and some of the surface
+is only required for fitting plots in `eigenstrapping.fit` and some of the surface
 plotting functions. Future improvements will reduce the number of dependencies
 needed.
 
@@ -81,10 +81,9 @@ Optional dependencies
 In order to speed up calculation of eigenmodes, you can utilize ``scikit-sparse`` libraries
 to use Cholesky decomposition rather than LU decomposition. If these libraries are already
 installed, the functions in ``eigenstrapping`` will automagically recognize this and use these
-libraries, without the user specifying (as default behavior, this can be turned off. Refer to
-:ref:`Deriving eigenmodes <usage_geometry_eigenmodes>`).
+libraries, without the user specifying (as default behavior, this can be turned off.
 
-You can install these libraries by following the instructions `here <https://github.com/scikit-sparse/scikit-sparse>`_.
+You can install these libraries by following the instructions: `scikit-sparse <https://github.com/scikit-sparse/scikit-sparse>`_.
 
 Citing
 ------
@@ -103,7 +102,7 @@ Please also cite the papers for the method that we use to calculate eigenmodes o
 
 * BrainPrint: a discriminative characterization of brain morphology. Wachinger C, Golland P, Kremen W, Fischl B, Reuter M. Neuroimage. 2015;109:232-48. `<http://dx.doi.org/10.1016/j.neuroimage.2015.01.032>`_ `<http://www.ncbi.nlm.nih.gov/pubmed/25613439>`_
 
-And if you use the Heat Kernel for Geodesics method in :func:`eigenstrapping.geometry.geodesic_distmat`:
+And if you use the Heat Kernel for Geodesics method in `eigenstrapping.geometry.geodesic_distmat`:
 
 * Crane, K., Weischedel, C., & Wardetzky, M. (2013). Geodesics in heat: A new approach to computing distance based on heat flow. ACM Transactions on Graphics (TOG), 32(5), 1-11. `<https://arxiv.org/pdf/1204.6216>`_
 

eigenstrapping/fit.py

@@ -19,7 +19,7 @@
                'distribution', 'shuffle', 'gen_rotations',
                'add_res']
 
-var_args = ['ns', 'pv', 'nh', 'knn', 'pv', 'nh', 'knn', 'n_jobs', 'seed']
+var_args = ['ns', 'pv', 'nh', 'knn', 'pv', 'nh', 'knn']
 
 def surface_fit(x, D=None, index=None, nsurrs=10, num_modes=100, return_data=False, 
                 extra_diags=False, surrs=None, **params):
@@ -93,7 +93,8 @@ def surface_fit(x, D=None, index=None, nsurrs=10, num_modes=100, return_data=Fal
         # surrogates
         surrs = eigen(n=nsurrs)
     else:
-        surrs = dataio(surrs)
+        surrs = surrs
+        nsurrs = len(surrs)
 
     # plot variogram
     # Instantiate surrogate map generator
@@ -175,7 +176,7 @@ def surface_fit(x, D=None, index=None, nsurrs=10, num_modes=100, return_data=Fal
     # Pairwise correlation plot
     ax = axes[1]
     correlations = efficient_pearsonr(x, surrs.T, nan_policy='omit')[0]
-    ax.hist(correlations, bins=np.linspace(-1, 1, num=nsurrs//5), color='#377eb8', alpha=0.7)
+    ax.hist(correlations, bins=np.linspace(-1, 1, num=nsurrs//10), color='#377eb8', alpha=0.7)
     ax.axvline(x=np.mean(correlations), color='r', linestyle='dashed', linewidth=3)
     ax.spines['top'].set_visible(False)
     ax.spines['right'].set_visible(False)

eigenstrapping/spatial.py

@@ -1,7 +1,8 @@
 # -*- coding: utf-8 -*-
 """
 Functions for calculating and manipulating spatial autocorrelation
-Author: Ross Markello, netneurolab
+Some functions: Ross Markello, netneurolab
+the rest: Nikitas Koussis, SNGNewy
 """
 
 import os
@@ -15,11 +16,16 @@
 from scipy.spatial import cKDTree
 from scipy.optimize import minimize
 from scipy.special import erfc
+from sklearn.utils.validation import check_random_state
 
 from netneurotools.datasets import make_correlated_xy
 from netneurotools.freesurfer import check_fs_subjid
 from netneurotools.utils import run
-from neuromaps import transforms, datasets
+from neuromaps import transforms, datasets, images
+
+from parspin import utils as putils
+
+from joblib import Parallel, delayed, dump, load
 
 from brainspace import mesh
 
@@ -718,20 +724,16 @@ def create_surface_grf(noise=None, alpha=3.0, normalize=True, seed=None,
     """
 
     affine = np.eye(4) * 2
-    affine[:, -1] = [-90, -90, -72, 1]
+    affine[:, -1] = [90, -126, -72, 1]
+    affine[0, 0] *= -1
 
     gfield = gaussian_random_field(91, 109, 91, noise=noise, alpha=alpha,
                                    normalize=normalize, seed=seed)
     fn = make_tmpname(suffix='.nii.gz')
-    nib.save(nib.nifti1.Nifti1Image(gfield, affine), fn)
+    nib.save(nib.nifti1.Nifti1Image(gfield.astype(np.float32), affine), fn)
 
-    data = np.zeros((20484,))
-    for n, hemi in enumerate(('lh', 'rh')):
-        outname = make_tmpname(suffix='.mgh')
-        run(VOL2SURF.format(fn, outname, hemi), quiet=True)
-        sl = slice(len(data) // 2 * n, len(data) // 2 * (n + 1))
-        data[sl] = nib.load(outname).get_fdata().squeeze()
-        os.remove(outname)
+    #data = np.zeros((20484,))
+    data = images.load_data(transforms.mni152_to_fsaverage(fn, fsavg_density='10k', method='nearest'))
 
     os.remove(fn)
 
@@ -740,6 +742,57 @@ def create_surface_grf(noise=None, alpha=3.0, normalize=True, seed=None,
 
     return data
 
+def calc_moran(dist, nulls, local=False, normalize=False, n_jobs=1, return_mask=False):
+    """
+    Calculates Moran's I for every column of `nulls`
+
+    Parameters
+    ----------
+    dist : (N, N) array_like
+        Full distance matrix (inter-hemispheric distance should be np.inf)
+    nulls : (N, P) array_like
+        Null brain maps for which to compute Moran's I
+    n_jobs : int, optional
+        Number of parallel workers to use for calculating Moran's I. Default: 1
+
+    Returns
+    -------
+    moran : (P,) np.ndarray
+        Moran's I for `P` null maps
+    """
+
+    def _moran(dist, sim, medmask, local=False, normalize=False, invert_dist=False):
+        mask = np.logical_and(medmask, np.logical_not(np.isnan(sim)))
+        return morans_i(dist[np.ix_(mask, mask)], sim[mask], local=local,
+                                normalize=normalize, invert_dist=invert_dist)
+
+    # do some pre-calculation on our distance matrix to reduce computation time
+    with np.errstate(divide='ignore', invalid='ignore'):
+        dist = 1 / dist
+        np.fill_diagonal(dist, 0)
+        dist /= dist.sum(axis=-1, keepdims=True)
+    # NaNs in the `dist` array are the "original" medial wall; mask these
+    medmask = np.logical_not(np.isnan(dist[:, 0]))
+
+    # calculate moran's I, masking out NaN values for each null (i.e., the
+    # rotated medial wall)
+    fn = dump(dist, make_tmpname('.mmap'))[0]
+    dist = load(fn, mmap_mode='r')
+    moran = np.array(
+        Parallel(n_jobs=n_jobs)(
+            delayed(_moran)(dist, nulls[:, n], medmask, local, normalize)
+            for n in putils.trange(nulls.shape[-1], desc="Running Moran's I")
+        )
+    )
+
+    mask = np.logical_and(medmask, np.logical_not(np.isnan(nulls[:, 0])))
+    Path(fn).unlink()
+
+    if return_mask:
+        return moran, mask
+
+    return moran
+
 def create_decimated_grf(surface, hemi='L', noise=None, alpha=3.0, normalize=True, seed=None,
                          medial=None):
     """
@@ -1002,6 +1055,51 @@ def matching_multinorm_grfs(corr, tol=0.005, *, alpha=3.0, normalize=True,
 
     return _mod_medial(xs, remove=False), _mod_medial(ys, remove=False)
 
+def grfs(alpha=3.0, normalize=True, seed=None, debug=False):
+    """
+    Generates two surface GRFs (fsaverage5)
+
+    Starts by generating two random variables from a multivariate normal
+    distribution, adds spatial autocorrelation with specified `alpha`, and 
+    projects to the surface.
+
+    Parameters
+    ----------
+    alpha : float (positive), optional
+        Exponent of the power-law distribution. Only used if `use_gstools` is
+        set to False. Default: 3.0
+    normalize : bool, optional
+        Whether to normalize the returned field to unit variance. Default: True
+    seed : None, int, default_rng, optional
+        Random state to seed GRF generation. Default: None
+    debug : bool, optional
+        Whether to print debug info
+
+    Return
+    ------
+    x, y : (20484,) np.ndarray
+        Generated surface GRFs
+    """
+
+    rs = np.random.default_rng(seed)
+
+    if alpha > 0:
+        # smooth correlated noise vectors + project to surface
+        xs = create_surface_grf(alpha=alpha, normalize=normalize)
+        ys = create_surface_grf(alpha=alpha, normalize=normalize)
+    else:
+        xs = rs.normal(size=20484)
+        ys = rs.normal(size=20484)
+
+        # remove medial wall to ensure data are still sufficiently correlated.
+        # this is important for parcellations that will ignore the medial wall
+        #xs, ys = _mod_medial(xs, remove=True), _mod_medial(ys, remove=True)
+
+    if normalize:
+        xs, ys = sstats.zscore(xs), sstats.zscore(ys)
+
+    return xs, ys #_mod_medial(xs, remove=False), _mod_medial(ys, remove=False)
+
 def matching_multinorm_exgrfs(corr, tol=0.005, alpha=3.0, params=[0., 1., 1.], normalize=True, seed=None, debug=False):
     """
     Generates two surface ex-GRFs (fsaverage5) that correlate at r = `corr`

eigenstrapping/stats.py

@@ -90,6 +90,14 @@ def msle(y_true, y_pred):
     """Compute the Mean Squared Logarithmic Error."""
     return np.mean((np.log1p(y_true) - np.log1p(y_pred)) ** 2)
 
+def fdr_bh(pvals):
+    """Compute the Benjamini-Hochberg procedure for false discovery rate correction"""
+    ranked_pvals = sstats.rankdata(pvals)
+    fdr = pvals * len(pvals) / ranked_pvals
+    fdr[fdr > 1] = 1
+    return fdr
+
+
 def joint_differential_entropy(x, y, b=0.5):
     """ Compute the joint differential entropy of x and y """
     data = np.vstack([x, y]).T

pyproject.toml

@@ -11,9 +11,8 @@ install_requires = [
 
 [project]
 name = "eigenstrapping"
-version = "0.0.2"
+version = "0.1"
 description = "For generating surrogate brain maps with spatial autocorrelation using geometric eigenmodes."
-readme = "README.rst"
 authors = [
     {name = "Nikitas C. Koussis", email = "[email protected]"},
 ]
@@ -30,7 +29,8 @@ dependencies = [
   "scipy>=0.11",
   "seaborn",
   "netneurotools>=0.2.4",
-  "tqdm"
+  "tqdm",
+  "parspin"
 ]
 maintainers = [
   {name = "Systems Neuroscience Group Newcastle", email = "[email protected]"},

requirements.txt

@@ -10,4 +10,5 @@ pandas
 scipy>=0.11
 seaborn
 netneurotools>=0.2.4
-tqdm
+tqdm
+parspin

setup.py

@@ -6,7 +6,7 @@
 
 # run the setup
 setup(name='eigenstrapping',
-  version='0.0.2',
+  version='0.1',
   description="For generating surrogate brain maps with spatial autocorrelation using geometric eigenmodes.",
   author='Nikitas C. Koussis, Systems Neuroscience Group Newcastle',
   author_email='[email protected]',