diff --git a/.circleci/config.yml b/.circleci/config.yml
index bb21645..bbcaac4 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -95,6 +95,74 @@ jobs:
           paths:
             - src/coverage/.coverage.py310
 
+  unittest_311:
+    docker:
+      - image: continuumio/miniconda3
+    working_directory: /tmp/src/mapca
+    steps:
+      - checkout
+      - restore_cache:
+          key: conda-py311-v1-{{ checksum "pyproject.toml" }}
+      - run:
+          name: Generate environment
+          command: |
+            apt-get update
+            apt-get install -yqq make
+            if [ ! -d /opt/conda/envs/mapca_py311 ]; then
+              conda create -yq -n mapca_py311 python=3.11
+              source activate mapca_py311
+              pip install .[tests]
+            fi
+      - run:
+          name: Running unit tests
+          command: |
+            source activate mapca_py311
+            py.test --ignore mapca/tests/test_integration.py --cov-append --cov-report term-missing --cov=mapca mapca/
+            mkdir /tmp/src/coverage
+            mv /tmp/src/mapca/.coverage /tmp/src/coverage/.coverage.py311
+      - save_cache:
+          key: conda-py311-v1-{{ checksum "pyproject.toml" }}
+          paths:
+            - /opt/conda/envs/mapca_py311
+      - persist_to_workspace:
+          root: /tmp
+          paths:
+            - src/coverage/.coverage.py311
+
+  unittest_312:
+    docker:
+      - image: continuumio/miniconda3
+    working_directory: /tmp/src/mapca
+    steps:
+      - checkout
+      - restore_cache:
+          key: conda-py312-v1-{{ checksum "pyproject.toml" }}
+      - run:
+          name: Generate environment
+          command: |
+            apt-get update
+            apt-get install -yqq make
+            if [ ! -d /opt/conda/envs/mapca_py312 ]; then
+              conda create -yq -n mapca_py312 python=3.12
+              source activate mapca_py312
+              pip install .[tests]
+            fi
+      - run:
+          name: Running unit tests
+          command: |
+            source activate mapca_py312
+            py.test --ignore mapca/tests/test_integration.py --cov-append --cov-report term-missing --cov=mapca mapca/
+            mkdir /tmp/src/coverage
+            mv /tmp/src/mapca/.coverage /tmp/src/coverage/.coverage.py312
+      - save_cache:
+          key: conda-py312-v1-{{ checksum "pyproject.toml" }}
+          paths:
+            - /opt/conda/envs/mapca_py312
+      - persist_to_workspace:
+          root: /tmp
+          paths:
+            - src/coverage/.coverage.py312
+
   unittest_38:
     docker:
       - image: continuumio/miniconda3
@@ -202,12 +270,13 @@ jobs:
             twine upload dist/*
 
 workflows:
-  version: 2.1
   build_test:
     jobs:
       - makeenv_38
       - unittest_39
       - unittest_310
+      - unittest_311
+      - unittest_312
       - unittest_38:
           requires:
             - makeenv_38
@@ -222,6 +291,8 @@ workflows:
             - unittest_38
             - unittest_39
             - unittest_310
+            - unittest_311
+            - unittest_312
       - deploy:
           requires:
             - merge_coverage
diff --git a/.gitignore b/.gitignore
index 8648d5e..14fd6fb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,6 @@
+# version file autogenerated with pip install -e .'[all]'
+_version.py
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]
diff --git a/mapca/__init__.py b/mapca/__init__.py
index 71868c6..a371305 100644
--- a/mapca/__init__.py
+++ b/mapca/__init__.py
@@ -1,13 +1,9 @@
-"""
-mapca: A Python implementation of the moving average principal components analysis methods from
-GIFT.
-"""
+"""mapca: Python implementation of the moving average principal components analysis from GIFT."""
 
 import warnings
 
 from mapca.__about__ import __version__
-
-from .mapca import MovingAveragePCA, ma_pca
+from mapca.mapca import MovingAveragePCA, ma_pca
 
 # cmp is not used, so ignore nipype-generated warnings
 warnings.filterwarnings("ignore", r"cmp not installed")
diff --git a/mapca/mapca.py b/mapca/mapca.py
index df7cbe4..a2fb43b 100644
--- a/mapca/mapca.py
+++ b/mapca/mapca.py
@@ -27,7 +27,7 @@
 from sklearn.decomposition import PCA
 from sklearn.preprocessing import StandardScaler
 
-from . import utils
+from mapca import utils
 
 LGR = logging.getLogger(__name__)
 
@@ -62,7 +62,7 @@ class MovingAveragePCA:
     explained_variance_ : :obj:`numpy.ndarray`, shape (n_components,)
         The amount of variance explained by each of the selected components.
 
-        Equal to n_components largest eigenvalues of the covariance matrix of X.
+        Equal to n_components largest eigenvalues of the covariance matrix of x.
     explained_variance_ratio_ : :obj:`numpy.ndarray`, shape (n_components,)
         Percentage of variance explained by each of the selected components.
 
@@ -75,7 +75,7 @@ class MovingAveragePCA:
     mean_ : :obj:`numpy.ndarray`, shape (n_features,)
         Per-feature empirical mean, estimated from the training set.
 
-        Equal to X.mean(axis=0).
+        Equal to x.mean(axis=0).
     n_components_ : int
         The estimated number of components.
         When n_components is set to ‘mle’ or a number between 0 and 1
@@ -144,31 +144,31 @@ def _fit(self, img, mask, subsample_depth=None):
         mask = mask.get_fdata()
 
         [n_x, n_y, n_z, n_timepoints] = data.shape
-        data_nib_V = np.reshape(data, (n_x * n_y * n_z, n_timepoints), order="F")
+        data_nib_v = np.reshape(data, (n_x * n_y * n_z, n_timepoints), order="F")
         mask_vec = np.reshape(mask, n_x * n_y * n_z, order="F")
-        X = data_nib_V[mask_vec == 1, :]
+        x = data_nib_v[mask_vec == 1, :]
 
         n_samples = np.sum(mask_vec)
 
         self.scaler_ = StandardScaler(with_mean=True, with_std=True)
         if self.normalize:
             # TODO: determine if tedana is already normalizing before this
-            X = self.scaler_.fit_transform(X.T).T  # This was X_sc
-            # X = ((X.T - X.T.mean(axis=0)) / X.T.std(axis=0)).T
+            x = self.scaler_.fit_transform(x.T).T  # This was x_sc
+            # x = ((x.T - x.T.mean(axis=0)) / x.T.std(axis=0)).T
 
         LGR.info("Performing SVD on original data...")
-        V, eigenvalues = utils._icatb_svd(X, n_timepoints)
+        v, eigenvalues = utils._icatb_svd(x, n_timepoints)
         LGR.info("SVD done on original data")
 
         # Reordering of values
         eigenvalues = eigenvalues[::-1]
-        dataN = np.dot(X, V[:, ::-1])
-        # Potentially the small differences come from the different signs on V
+        data_n = np.dot(x, v[:, ::-1])
+        # Potentially the small differences come from the different signs on v
 
         # Using 12 gaussian components from middle, top and bottom gaussian
         # components to determine the subsampling depth.
         # Final subsampling depth is determined using median
-        kurt = kurtosis(dataN, axis=0, fisher=True)
+        kurt = kurtosis(data_n, axis=0, fisher=True)
         kurt[kurt < 0] = 0
         kurt = np.expand_dims(kurt, 1)
 
@@ -180,25 +180,25 @@ def _fit(self, img, mask, subsample_depth=None):
         ]  # NOTE: make sure np.where is giving us just one tuple
         idx = np.array(idx_gauss[:]).T
         dfs = np.sum(eigenvalues > np.finfo(float).eps)  # degrees of freedom
-        minTp = 12
+        min_tp = 12
 
-        if len(idx) >= minTp:
+        if len(idx) >= min_tp:
             middle = int(np.round(len(idx) / 2))
             idx = np.hstack([idx[0:4], idx[middle - 1 : middle + 3], idx[-4:]])
         else:
-            minTp = np.min([minTp, dfs])
-            idx = np.arange(dfs - minTp, dfs)
+            min_tp = np.min([min_tp, dfs])
+            idx = np.arange(dfs - min_tp, dfs)
 
         idx = np.unique(idx)
 
         # Estimate the subsampling depth for effectively i.i.d. samples
         LGR.info("Estimating the subsampling depth for effective i.i.d samples...")
-        mask_ND = np.reshape(mask_vec, (n_x, n_y, n_z), order="F")
+        mask_nd = np.reshape(mask_vec, (n_x, n_y, n_z), order="F")
         sub_depth = len(idx)
         sub_iid_sp = np.zeros((sub_depth,))
         for i in range(sub_depth):
             x_single = np.zeros(n_x * n_y * n_z)
-            x_single[mask_vec == 1] = dataN[:, idx[i]]
+            x_single[mask_vec == 1] = data_n[:, idx[i]]
             x_single = np.reshape(x_single, (n_x, n_y, n_z), order="F")
             sub_iid_sp[i] = utils._est_indp_sp(x_single)[0] + 1
             if i > 6:
@@ -223,7 +223,7 @@ def _fit(self, img, mask, subsample_depth=None):
             )
             sub_iid_sp_median = int(np.floor(np.power(n_samples / n_timepoints, 1 / dim_n)))
 
-        LGR.info("Estimated subsampling depth for effective i.i.d samples: %d" % sub_iid_sp_median)
+        LGR.info(f"Estimated subsampling depth for effective i.i.d samples: {sub_iid_sp_median}")
 
         # Always save the calculated IID subsample value, but, if there is a user provide value,
         # assign that to sub_iid_sp_median and use that instead
@@ -250,19 +250,19 @@ def _fit(self, img, mask, subsample_depth=None):
                 # probably good to at least include an upper bound.
                 raise ValueError(
                     "subsample_depth must be an integer > 1 and will retain >100 "
-                    "samples after subsampling. It is %d" % subsample_depth
+                    f"samples after subsampling. It is {subsample_depth}"
                 )
 
-        N = np.round(n_samples / np.power(sub_iid_sp_median, dim_n))
+        n = np.round(n_samples / np.power(sub_iid_sp_median, dim_n))
 
         if sub_iid_sp_median != 1:
-            mask_s = utils._subsampling(mask_ND, sub_iid_sp_median)
+            mask_s = utils._subsampling(mask_nd, sub_iid_sp_median)
             mask_s_1d = np.reshape(mask_s, np.prod(mask_s.shape), order="F")
             dat = np.zeros((int(np.sum(mask_s_1d)), n_timepoints))
             LGR.info("Generating subsampled i.i.d. data...")
             for i_vol in range(n_timepoints):
                 x_single = np.zeros(n_x * n_y * n_z)
-                x_single[mask_vec == 1] = X[:, i_vol]
+                x_single[mask_vec == 1] = x[:, i_vol]
                 x_single = np.reshape(x_single, (n_x, n_y, n_z), order="F")
                 dat0 = utils._subsampling(x_single, sub_iid_sp_median)
                 dat0 = np.reshape(dat0, np.prod(dat0.shape), order="F")
@@ -274,15 +274,15 @@ def _fit(self, img, mask, subsample_depth=None):
 
             # (completed)
             LGR.info("Performing SVD on subsampled i.i.d. data...")
-            V, eigenvalues = utils._icatb_svd(dat, n_timepoints)
+            v, eigenvalues = utils._icatb_svd(dat, n_timepoints)
             LGR.info("SVD done on subsampled i.i.d. data")
             eigenvalues = eigenvalues[::-1]
 
-        LGR.info("Effective number of i.i.d. samples %d from %d total voxels" % (N, n_samples))
+        LGR.info(f"Effective number of i.i.d. samples {n} from {n_samples} total voxels")
 
         # Make eigen spectrum adjustment
         LGR.info("Perform eigen spectrum adjustment ...")
-        eigenvalues = utils._eigensp_adj(eigenvalues, N, eigenvalues.shape[0])
+        eigenvalues = utils._eigensp_adj(eigenvalues, n, eigenvalues.shape[0])
         # (completed)
         if np.sum(np.imag(eigenvalues)):
             raise ValueError("Invalid eigenvalue found for the subsampled data.")
@@ -301,12 +301,12 @@ def _fit(self, img, mask, subsample_depth=None):
         mdl = np.zeros(p - 1)
 
         for k_idx, k in enumerate(np.arange(1, p)):
-            LH = np.log(np.prod(np.power(eigenvalues[k:], 1 / (p - k))) / np.mean(eigenvalues[k:]))
-            mlh = 0.5 * N * (p - k) * LH
+            lh = np.log(np.prod(np.power(eigenvalues[k:], 1 / (p - k))) / np.mean(eigenvalues[k:]))
+            mlh = 0.5 * n * (p - k) * lh
             df = 1 + 0.5 * k * (2 * p - k + 1)
             aic[k_idx] = (-2 * mlh) + (2 * df)
             kic[k_idx] = (-2 * mlh) + (3 * df)
-            mdl[k_idx] = -mlh + (0.5 * df * np.log(N))
+            mdl[k_idx] = -mlh + (0.5 * df * np.log(n))
 
         itc = np.row_stack([aic, kic, mdl])
 
@@ -345,7 +345,7 @@ def _fit(self, img, mask, subsample_depth=None):
 
         # PCA with all possible components (the estimated selection is made after)
         ppca = PCA(n_components=None, svd_solver="full", copy=False, whiten=False)
-        ppca.fit(X)
+        ppca.fit(x)
 
         # Get cumulative explained variance as components are added
         cumsum_varexp = np.cumsum(ppca.explained_variance_ratio_)
@@ -356,7 +356,7 @@ def _fit(self, img, mask, subsample_depth=None):
         # Calculate number of components for 95% varexp
         n_comp_varexp_95 = np.where(cumsum_varexp >= 0.95)[0][0] + 1
 
-        LGR.info("Estimated number of components is %d" % n_components)
+        LGR.info(f"Estimated number of components is {n_components}")
 
         # Save results of each criterion into dictionaries
         self.aic_ = {
@@ -391,7 +391,7 @@ def _fit(self, img, mask, subsample_depth=None):
             "calculated_IID_subsample_mean": sub_iid_sp_mean,
             "IID_subsample_input": sub_iid_sp,
             "used_IID_subsample_depth": sub_iid_sp_median,
-            "effective_num_IID_samples": N,
+            "effective_num_IID_samples": n,
             "total_num_samples": n_samples,
         }
 
@@ -407,7 +407,7 @@ def _fit(self, img, mask, subsample_depth=None):
         # Commenting out noise variance as it depends on the covariance of the estimation
         # self.noise_variance_ = ppca.noise_variance_
         component_maps = np.dot(
-            np.dot(X, self.components_.T), np.diag(1.0 / self.explained_variance_)
+            np.dot(x, self.components_.T), np.diag(1.0 / self.explained_variance_)
         )
         component_maps_3d = np.zeros((n_x * n_y * n_z, n_components))
         component_maps_3d[mask_vec == 1, :] = component_maps
@@ -416,7 +416,7 @@ def _fit(self, img, mask, subsample_depth=None):
         self.u_nii_ = nib.Nifti1Image(component_maps_3d, img.affine, img.header)
 
     def fit(self, img, mask, subsample_depth=None):
-        """Fit the model with X.
+        """Fit the model with x.
 
         Parameters
         ----------
@@ -442,7 +442,7 @@ def fit(self, img, mask, subsample_depth=None):
         return self
 
     def fit_transform(self, img, mask, subsample_depth=None):
-        """Fit the model with X and apply the dimensionality reduction on X.
+        """Fit the model with x and apply the dimensionality reduction on x.
 
         Parameters
         ----------
@@ -458,7 +458,7 @@ def fit_transform(self, img, mask, subsample_depth=None):
 
         Returns
         -------
-        X_new : 4D niimg_like
+        x_new : 4D niimg_like
             Component weight maps.
 
         Notes
@@ -477,34 +477,25 @@ def fit_transform(self, img, mask, subsample_depth=None):
         dataset.
         """
         self._fit(img, mask, subsample_depth=subsample_depth)
-        return self.transform(img)
+        return self.transform()
 
-    def transform(self, img):
-        """Apply dimensionality reduction to X.
-
-        X is projected on the first principal components previously extracted from a training set.
-
-        Parameters
-        ----------
-        img : 4D niimg_like
-            Data on which to apply PCA.
+    def transform(self):
+        """Return x after dimensionality reduction.
 
         Returns
         -------
-        X_new : array-like, shape (n_samples, n_components)
+        x_new : array-like, shape (n_samples, n_components)
 
         Notes
         -----
         This is different from scikit-learn's approach, which ignores explained variance.
         """
-        # X = self.scaler_.fit_transform(X.T).T
-        # X_new = np.dot(np.dot(X, self.components_.T), np.diag(1.0 / self.explained_variance_))
         return self.u_nii_
 
     def inverse_transform(self, img, mask):
         """Transform data back to its original space.
 
-        In other words, return an input X_original whose transform would be X.
+        In other words, return an input x_original whose transform would be x.
 
         Parameters
         ----------
@@ -528,17 +519,17 @@ def inverse_transform(self, img, mask):
         mask = mask.get_fdata()
 
         [n_x, n_y, n_z, n_components] = data.shape
-        data_nib_V = np.reshape(data, (n_x * n_y * n_z, n_components), order="F")
+        data_nib_v = np.reshape(data, (n_x * n_y * n_z, n_components), order="F")
         mask_vec = np.reshape(mask, n_x * n_y * n_z, order="F")
-        X = data_nib_V[mask_vec == 1, :]
+        x = data_nib_v[mask_vec == 1, :]
 
-        X_orig = np.dot(np.dot(X, np.diag(self.explained_variance_)), self.components_)
+        x_orig = np.dot(np.dot(x, np.diag(self.explained_variance_)), self.components_)
         if self.normalize:
-            X_orig = self.scaler_.inverse_transform(X_orig.T).T
+            x_orig = self.scaler_.inverse_transform(x_orig.T).T
 
-        n_t = X_orig.shape[1]
+        n_t = x_orig.shape[1]
         out_data = np.zeros((n_x * n_y * n_z, n_t))
-        out_data[mask_vec == 1, :] = X_orig
+        out_data[mask_vec == 1, :] = x_orig
         out_data = np.reshape(out_data, (n_x, n_y, n_z, n_t), order="F")
         img_orig = nib.Nifti1Image(out_data, img.affine, img.header)
         return img_orig
diff --git a/mapca/tests/conftest.py b/mapca/tests/conftest.py
index facdffa..8020020 100644
--- a/mapca/tests/conftest.py
+++ b/mapca/tests/conftest.py
@@ -1,3 +1,5 @@
+"""Tests for mapca configuration."""
+
 import os
 from urllib.request import urlretrieve
 
@@ -6,7 +8,7 @@
 
 def fetch_file(osf_id, path, filename):
     """
-    Fetches file located on OSF and downloads to `path`/`filename`1
+    Fetch file located on OSF and downloads to `path`/`filename`1.
 
     Parameters
     ----------
@@ -25,7 +27,7 @@ def fetch_file(osf_id, path, filename):
     full_path : str
         Full path to downloaded `filename`
     """
-    url = "https://osf.io/{}/download".format(osf_id)
+    url = f"https://osf.io/{osf_id}/download"
     full_path = os.path.join(path, filename)
     if not os.path.isfile(full_path):
         urlretrieve(url, full_path)
@@ -34,35 +36,41 @@ def fetch_file(osf_id, path, filename):
 
 @pytest.fixture(scope="session")
 def testpath(tmp_path_factory):
-    """Test path that will be used to download all files"""
+    """Test path that will be used to download all files."""
     return tmp_path_factory.getbasetemp()
 
 
 @pytest.fixture
 def test_img(testpath):
+    """Fetch data file."""
     return fetch_file("jw43h", testpath, "data.nii.gz")
 
 
 @pytest.fixture
 def test_mask(testpath):
+    """Fetch mask file."""
     return fetch_file("9u2m5", testpath, "mask.nii.gz")
 
 
 @pytest.fixture
 def test_ts(testpath):
+    """Fetch comp_ts file."""
     return fetch_file("gz2hb", testpath, "comp_ts.npy")
 
 
 @pytest.fixture
 def test_varex(testpath):
+    """Fetch varex file."""
     return fetch_file("7xj5k", testpath, "varex.npy")
 
 
 @pytest.fixture
 def test_varex_norm(testpath):
+    """Fetch varex_norm file."""
     return fetch_file("jrd9c", testpath, "varex_norm.npy")
 
 
 @pytest.fixture
 def test_weights(testpath):
+    """Fetch weights file."""
     return fetch_file("t94m8", testpath, "voxel_comp_weights.npy")
diff --git a/mapca/tests/test_integration.py b/mapca/tests/test_integration.py
index 0a9ca30..3235d7f 100644
--- a/mapca/tests/test_integration.py
+++ b/mapca/tests/test_integration.py
@@ -1,6 +1,4 @@
-"""
-Integration test for mapca.
-"""
+"""Integration tests for mapca."""
 
 import shutil
 from os.path import split
diff --git a/mapca/tests/test_mapca.py b/mapca/tests/test_mapca.py
index 1dc88ad..9213085 100644
--- a/mapca/tests/test_mapca.py
+++ b/mapca/tests/test_mapca.py
@@ -1,6 +1,4 @@
-"""
-Tests for mapca
-"""
+"""Tests for mapca."""
 
 import nibabel as nib
 import numpy as np
@@ -10,8 +8,7 @@
 
 
 def test_ma_pca():
-    """Check that ma_pca runs correctly with all three options"""
-
+    """Check that ma_pca runs correctly with all three options."""
     n_timepoints = 200
     n_voxels = 20
     n_vox_total = n_voxels**3
@@ -57,22 +54,21 @@ def test_ma_pca():
     assert v.shape[0] == n_timepoints
 
 
-def test_MovingAveragePCA():
+def test_moving_average_pca():
     """Check that MovingAveragePCA runs correctly with "aic" option."""
-
-    N_TIMEPOINTS = 200
-    N_VOXELS = 20  # number of voxels in each dimension
+    n_timepoints = 200
+    n_voxels = 20  # number of voxels in each dimension
 
     # Create fake data to test with
-    test_data = np.random.random((N_VOXELS, N_VOXELS, N_VOXELS, N_TIMEPOINTS))
-    time = np.linspace(0, 400, N_TIMEPOINTS)
+    test_data = np.random.random((n_voxels, n_voxels, n_voxels, n_timepoints))
+    time = np.linspace(0, 400, n_timepoints)
     freq = 1
     test_data = test_data + np.sin(2 * np.pi * freq * time)
     xform = np.eye(4) * 2
     test_img = nib.nifti1.Nifti1Image(test_data, xform)
 
     # Create mask
-    test_mask = np.zeros((N_VOXELS, N_VOXELS, N_VOXELS), dtype=int)
+    test_mask = np.zeros((n_voxels, n_voxels, n_voxels), dtype=int)
     test_mask[5:-5, 5:-5, 5:-5] = 1
     test_mask_img = nib.nifti1.Nifti1Image(test_mask, xform, dtype=np.int16)
     n_voxels_in_mask = np.sum(test_mask)
@@ -85,12 +81,12 @@ def test_MovingAveragePCA():
     assert pca.u_.shape[0] == n_voxels_in_mask
     assert pca.explained_variance_.shape[0] == 1
     assert pca.explained_variance_ratio_.shape[0] == 1
-    assert pca.components_.T.shape[0] == N_TIMEPOINTS
+    assert pca.components_.T.shape[0] == n_timepoints
 
     # Test other stuff
     pca2 = MovingAveragePCA(criterion="mdl", normalize=True)
     pca2.fit(test_img, test_mask_img)
-    u2 = pca2.transform(test_img)
+    u2 = pca2.transform()
     assert np.array_equal(u2.get_fdata(), u.get_fdata())
 
     test_data_est = pca2.inverse_transform(u2, test_mask_img)
diff --git a/mapca/tests/test_utils.py b/mapca/tests/test_utils.py
index bb8e354..ef02b1b 100644
--- a/mapca/tests/test_utils.py
+++ b/mapca/tests/test_utils.py
@@ -10,9 +10,7 @@
 
 
 def test_autocorr():
-    """
-    Unit test on _autocorr function
-    """
+    """Unit test on _autocorr function."""
     test_data = np.array([1, 2, 3, 4])
     test_result = np.array([30, 20, 11, 4])
     autocorr = _autocorr(test_data)
@@ -20,6 +18,7 @@ def test_autocorr():
 
 
 def test_parzen_win():
+    """Test parzen gives expected output."""
     test_npoints = 3
     test_result = np.array([0.07407407, 1, 0.07407407])
     win = parzen(test_npoints)
@@ -31,9 +30,7 @@ def test_parzen_win():
 
 
 def test_ent_rate_sp():
-    """
-    Check that ent_rate_sp runs correctly, i.e. returns a float
-    """
+    """Check that ent_rate_sp runs correctly, i.e. returns a float."""
     test_data = np.random.rand(200, 10, 10)
     ent_rate = ent_rate_sp(test_data, 0)
     ent_rate = ent_rate_sp(test_data, 1)
@@ -55,9 +52,7 @@ def test_ent_rate_sp():
 
 
 def test_subsampling():
-    """
-    Unit test for subsampling function
-    """
+    """Unit test for subsampling function."""
     # 1D input
     test_data = np.array([1])
     sub_data = _subsampling(test_data, sub_depth=2)
@@ -76,18 +71,14 @@ def test_subsampling():
 
 
 def test_icatb_svd():
-    """
-    Unit test for icatb_svd function.
-    """
+    """Unit test for icatb_svd function."""
     test_data = np.diag(np.random.rand(5))
-    V, Lambda = _icatb_svd(test_data)
-    assert np.allclose(np.sum(V, axis=0), np.ones((5,)))
+    v, lambda_var = _icatb_svd(test_data)
+    assert np.allclose(np.sum(v, axis=0), np.ones((5,)))
 
 
 def test_eigensp_adj():
-    """
-    Unit test for eigensp_adj function
-    """
+    """Unit test for eigensp_adj function."""
     test_eigen = np.array([0.9, 0.5, 0.2, 0.1, 0])
     n_effective = 2
     test_result = np.array([0.13508894, 0.11653465, 0.06727316, 0.05211424, 0.0])
@@ -96,7 +87,7 @@ def test_eigensp_adj():
 
 
 def test_kurtosis():
-    # Generate data
+    """Generate data."""
     test_data = np.array([[-10, 2, 500, 0, -0.4], [-4, -200, -40, 0.1, 90]]).T
 
     # Run scipy function
diff --git a/mapca/utils.py b/mapca/utils.py
index b682447..a2ebc17 100644
--- a/mapca/utils.py
+++ b/mapca/utils.py
@@ -1,6 +1,4 @@
-"""
-PCA based on Moving Average (stationary Gaussian) process
-"""
+"""PCA based on Moving Average (stationary Gaussian) process."""
 
 import logging
 
@@ -16,7 +14,7 @@
 
 def _autocorr(data):
     """
-    Calculates the auto correlation of a given array.
+    Calculate the auto correlation of a given array.
 
     Parameters
     ----------
@@ -35,8 +33,9 @@ def _autocorr(data):
 
 def ent_rate_sp(data, sm_window):
     """
-    Calculate the entropy rate of a stationary Gaussian random process using
-    spectrum estimation with smoothing window.
+    Calculate the entropy rate of a stationary Gaussian random process.
+
+    Uses spectrum estimation with smoothing window.
 
     Parameters
     ----------
@@ -61,7 +60,6 @@ def ent_rate_sp(data, sm_window):
       functional magnetic resonance imaging data.
       Human brain mapping, 28(11), pp.1251-1266.
     """
-
     dims = data.shape
 
     if data.ndim == 3 and min(dims) != 1:
@@ -94,30 +92,30 @@ def ent_rate_sp(data, sm_window):
     data_corr /= vcu
 
     if sm_window:
-        M = [int(i) for i in np.ceil(np.array(dims) / 10)]
+        m = [int(i) for i in np.ceil(np.array(dims) / 10)]
 
         # Get Parzen window for each spatial direction
         parzen_w_3 = np.zeros((2 * dims[2] - 1,))
-        parzen_w_3[(dims[2] - M[2] - 1) : (dims[2] + M[2])] = parzen(2 * M[2] + 1)
+        parzen_w_3[(dims[2] - m[2] - 1) : (dims[2] + m[2])] = parzen(2 * m[2] + 1)
 
         parzen_w_2 = np.zeros((2 * dims[1] - 1,))
-        parzen_w_2[(dims[1] - M[1] - 1) : (dims[1] + M[1])] = parzen(2 * M[1] + 1)
+        parzen_w_2[(dims[1] - m[1] - 1) : (dims[1] + m[1])] = parzen(2 * m[1] + 1)
 
         parzen_w_1 = np.zeros((2 * dims[0] - 1,))
-        parzen_w_1[(dims[0] - M[0] - 1) : (dims[0] + M[0])] = parzen(2 * M[0] + 1)
+        parzen_w_1[(dims[0] - m[0] - 1) : (dims[0] + m[0])] = parzen(2 * m[0] + 1)
 
         # Scale Parzen windows
-        parzen_window_2D = np.dot(parzen_w_1[np.newaxis, :].T, parzen_w_2[np.newaxis, :])
-        parzen_window_3D = np.zeros((2 * dims[0] - 1, 2 * dims[1] - 1, 2 * dims[2] - 1))
+        parzen_window_2d = np.dot(parzen_w_1[np.newaxis, :].T, parzen_w_2[np.newaxis, :])
+        parzen_window_3d = np.zeros((2 * dims[0] - 1, 2 * dims[1] - 1, 2 * dims[2] - 1))
         for m3 in range(dims[2] - 1):
-            parzen_window_3D[:, :, (dims[2] - 1) - m3] = np.dot(
-                parzen_window_2D, parzen_w_3[dims[2] - 1 - m3]
+            parzen_window_3d[:, :, (dims[2] - 1) - m3] = np.dot(
+                parzen_window_2d, parzen_w_3[dims[2] - 1 - m3]
             )
-            parzen_window_3D[:, :, (dims[2] - 1) + m3] = np.dot(
-                parzen_window_2D, parzen_w_3[dims[2] - 1 + m3]
+            parzen_window_3d[:, :, (dims[2] - 1) + m3] = np.dot(
+                parzen_window_2d, parzen_w_3[dims[2] - 1 + m3]
             )
         # Apply 3D Parzen Window
-        data_corr *= parzen_window_3D
+        data_corr *= parzen_window_3d
 
     data_fft = abs(fftshift(fftn(data_corr)))
     data_fft[data_fft < 1e-4] = 1e-4
@@ -132,8 +130,9 @@ def ent_rate_sp(data, sm_window):
 
 def _est_indp_sp(data):
     """
-    Estimate the effective number of independent samples based on the maximum
-    entropy rate principle of stationary random process.
+    Estimate the effective number of independent samples.
+
+    Bbased on the maximum entropy rate principle of stationary random process.
 
     Parameters
     ----------
@@ -152,7 +151,6 @@ def _est_indp_sp(data):
     This function estimates the effective number of independent samples by omitting
     the least significant components with the subsampling scheme (Li et al., 2007)
     """
-
     dims = data.shape
     n_iters_0 = None
 
@@ -173,8 +171,7 @@ def _est_indp_sp(data):
         raise ValueError("Ill conditioned data, can not estimate " "independent samples.")
     n_iters = n_iters_0
     LGR.debug(
-        "Estimated the entropy rate of the Gaussian component "
-        "with subsampling depth {}".format(j + 1)
+        f"Estimated the entropy rate of the Gaussian component with subsampling depth {j + 1}"
     )
 
     return n_iters, ent_rate
@@ -196,7 +193,6 @@ def _subsampling(data, sub_depth):
     out : ndarray
         Subsampled data
     """
-
     slices = [slice(None, None, sub_depth)] * data.ndim
     out = data[tuple(slices)]
     return out
@@ -204,7 +200,7 @@ def _subsampling(data, sub_depth):
 
 def _kurtn(data):
     """
-    Normalized kurtosis funtion so that for a Gaussian r.v. the kurtn(g) = 0.
+    Normalize kurtosis funtion so that for a Gaussian r.v. the kurtn(g) = 0.
 
     Parameters
     ----------
@@ -217,7 +213,6 @@ def _kurtn(data):
         The kurtosis of each vector in x along the second dimension. For
         tedana, this will be the kurtosis of each PCA component.
     """
-
     kurt = kurtosis(data, axis=0, fisher=True)
     kurt[kurt < 0] = 0
     kurt = np.expand_dims(kurt, 1)
@@ -227,8 +222,9 @@ def _kurtn(data):
 
 def _icatb_svd(data, n_comps=None):
     """
-    Run Singular Value Decomposition (SVD) on input data and extracts the
-    given number of components (n_comps).
+    Run Singular Value Decomposition (SVD) on input data.
+
+    Runs SVD and extracts the given number of components (n_comps).
 
     Parameters
     ----------
@@ -244,28 +240,27 @@ def _icatb_svd(data, n_comps=None):
     Lambda : float
         Eigenvalues
     """
-
     if not n_comps:
         n_comps = np.min((data.shape[0], data.shape[1]))
 
-    _, Lambda, vh = svd(data, full_matrices=False)
+    _, lambda_var, vh = svd(data, full_matrices=False)
 
     # Sort eigen vectors in Ascending order
-    V = vh.T
-    Lambda = Lambda / np.sqrt(data.shape[0] - 1)  # Whitening (sklearn)
-    inds = np.argsort(np.power(Lambda, 2))
-    Lambda = np.power(Lambda, 2)[inds]
-    V = V[:, inds]
-    sumAll = np.sum(Lambda)
+    v = vh.T
+    lambda_var = lambda_var / np.sqrt(data.shape[0] - 1)  # Whitening (sklearn)
+    inds = np.argsort(np.power(lambda_var, 2))
+    lambda_var = np.power(lambda_var, 2)[inds]
+    v = v[:, inds]
+    sum_all = np.sum(lambda_var)
 
     # Return only the extracted components
-    V = V[:, (V.shape[1] - n_comps) :]
-    Lambda = Lambda[Lambda.shape[0] - n_comps :]
-    sumUsed = np.sum(Lambda)
-    retained = (sumUsed / sumAll) * 100
-    LGR.debug("{ret}% of non-zero components retained".format(ret=retained))
+    v = v[:, (v.shape[1] - n_comps) :]
+    lambda_var = lambda_var[lambda_var.shape[0] - n_comps :]
+    sum_used = np.sum(lambda_var)
+    retained = (sum_used / sum_all) * 100
+    LGR.debug(f"{retained}% of non-zero components retained")
 
-    return V, Lambda
+    return v, lambda_var
 
 
 def _eigensp_adj(lam, n, p):
@@ -298,7 +293,6 @@ def _eigensp_adj(lam, n, p):
       functional magnetic resonance imaging data.
       Human brain mapping, 28(11), pp.1251-1266.
     """
-
     r = p / n
     bp = np.power((1 + np.sqrt(r)), 2)
     bm = np.power((1 - np.sqrt(r)), 2)
diff --git a/pyproject.toml b/pyproject.toml
index 03b5b70..0f32581 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -15,6 +15,9 @@ classifiers = [
     "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",
+
 ]
 license = {file = "LICENSE"}
 requires-python = ">=3.8"