original commit

denoise/config.yml

@@ -0,0 +1,88 @@
+regr_set:
+  - vc
+  - fm
+  - dm
+  - mp_pca
+  - ica
+  - bp
+
+defaults:
+  vc: 
+    structure:
+      - derivatives
+    name:
+      - task-rest_desc-censored_vols.txt
+    transform_method:
+      - vol_censor
+    ismap:
+      - False
+    options:
+      - null
+  fm:
+    structure:
+      - rawdata
+    name:
+      - task-rest_rec-denoisemapfold1.nii.gz
+      - task-rest_rec-denoisemapfold2.nii.gz
+    transform_method:
+      - zscoring
+    ismap:
+      - True
+    options:
+      - null
+  dm:
+    structure:
+      - rawdata
+    name:
+      - task-rest_rec-denoisemapdens.nii.gz
+    transform_method:
+      - zscoring
+    ismap:
+      - True
+    options:
+      - null
+  mp_pca:
+    structure:
+      - rawdata
+    name:
+      - task-rest_rec-denoisemapxrot.nii.gz 
+      - task-rest_rec-denoisemapxtransl.nii.gz
+      - task-rest_rec-denoisemapyrot.nii.gz
+      - task-rest_rec-denoisemapytransl.nii.gz
+      - task-rest_rec-denoisemapzrot.nii.gz
+      - task-rest_rec-denoisemapztransl.nii.gz
+      - task-rest_rec-denoisemapdxrot.nii.gz 
+      - task-rest_rec-denoisemapdxtransl.nii.gz
+      - task-rest_rec-denoisemapdyrot.nii.gz
+      - task-rest_rec-denoisemapdytransl.nii.gz
+      - task-rest_rec-denoisemapdzrot.nii.gz
+      - task-rest_rec-denoisemapdztransl.nii.gz
+    transform_method:
+      - mp_pca
+    ismap:
+      - True
+    options:
+      - 36_24_p99
+  ica:
+    structure:
+      - derivatives
+    name:
+      - task-rest_desc-wm_csf_timeseries
+    transform_method:
+      - zscoring
+    ismap:
+      - False
+    options:
+      - 6
+  bp:
+    structure:
+      - null
+    name:
+      - null
+    transform_method:
+      - bptf
+    ismap:
+      - False
+    options:
+      - h0067
+

denoise/dnx_wrapper.py

@@ -0,0 +1,178 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+    @authors: Slava Karolis (slava.karolis@kcl.ac.uk) & Denis Prokopenko (denis.prokopenko@kcl.ac.uk)
+    
+    This work is jointly copyrighted by University of Oxford and King's College London
+    Copyright 2024
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+        http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+"""
+from os import getcwd
+import os
+import sys
+
+sys.path.append(getcwd())
+from methods.dnx import DV, IV, Regress
+import argparse
+import yaml
+from ios.vol_io import save_vol
+
+
+def sort_out_options(regr_type, options_string):
+    if regr_type == "mp_pca":
+        options_split = options_string.split(
+            "_"
+        )  # eg '24_12_p99' - 24 predictors from each subset and cut at 99th percentile
+        options = []
+        for i in range(2):
+            nregs = int(options_split[i])
+            mod = (
+                nregs % 6
+            )  # ie. needs to conform to the fact that the number of mp regressors == 6
+            if mod == 0:
+                options.append(nregs)
+            else:
+                if i == 0:
+                    errmsg = "The value for the number of volume-to-volume mp regressors has to be devisible by 6 without a remainder"
+                else:
+                    errmsg = "The value for the number of slices-to-slices mp regressors has to be devisible by 6 without a remainder"
+                raise ValueError(errmsg)
+
+        if options_split[2][0] == "p":
+            options.append(
+                float(options_split[2][1:]) / 100
+            )  # convert to the value between 0 and 1 (proportion of variance accounted for by pca)
+        elif options_split[2][0] == "n":
+            options.append(int(options_split[2][1:]))
+        else:
+            raise "Wrong input for the mp-based regressor threshold, needs to be prepended with n (#components) or p (variance accounted)"
+    elif regr_type == "bp":
+        # options format needs to be in "l1_h01" format  == lowpass =.1 & highpass==.01
+        # or l1 or hp
+
+        options_split = options_string.split("_")
+
+        for i in options_split:
+            if i[0] == "l":
+                lp = list(i)
+                lp[0] = "."
+                lp = float("".join(lp))
+            else:
+                hp = list(i)
+                hp[0] = "."
+                hp = float("".join(hp))
+        options = []
+        if "lp" not in locals():
+            options.append(None)
+        else:
+            options.append(lp)
+        if "hp" not in locals():
+            options.append(None)
+        else:
+            options.append(hp)
+    else:
+        if options_string == "None":
+            options = [None]
+        elif options_string is None:
+            options = [options_string]
+
+        else:
+            options = [int(options_string)]
+    return options
+
+
+def wrapper(root_folder, sub, ses, config):
+
+    dv = DV(sub, ses, root_folder)
+    iv = IV(sub, ses, root_folder, config)
+    rgx = Regress()
+    dv.activate_dnx(iv, rgx)
+
+    return dv, iv, rgx
+
+
+def save_dnx_vols(dv, odir):
+    if os.path.exists(odir) == False:
+        os.mkdir(odir)
+
+    # residuals
+    save_vol(dv.resid, odir + "/data.nii.gz", dv.mask_path, dv.path)
+    # yhat
+    save_vol(dv.yhat, odir + "/yhat.nii.gz", dv.mask_path, dv.path)
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(
+        prog="ProgName",
+        description="What the program does",
+        epilog="Text at the bottom of help",
+    )
+
+    parser.add_argument(
+        "-c", "--config", help="Path to config file", default="./config.yml"
+    )
+    parser.add_argument("--sub", help="Subject ID", required=True)
+    parser.add_argument("--ses", help="Session ID", required=True)
+    parser.add_argument(
+        "--odir", help="directory for writing the output", required=True
+    )
+    parser.add_argument(
+        "--root_dir", help="root directory for the dHCP dataset", required=True
+    )
+    parser.add_argument(
+        "--regr_set",
+        nargs="+",
+        help="regressor set, using space as a separator",
+        required=False,
+    )
+    parser.add_argument(
+        "--ica-options",
+        help="user-defined number of components for ica-based regressors",
+        required=False,
+    )
+    parser.add_argument(
+        "--mp_pca-options",
+        help="user-defined options for mp-based regressors",
+        required=False,
+    )
+    parser.add_argument(
+        "--bp-options", help="user-defined temporal frequency bandpass", required=False
+    )
+
+    args = vars(parser.parse_args())
+    print(args)
+
+    with open(args["config"], "r") as f:
+        default_config = yaml.safe_load(f)
+
+    if args["regr_set"] is not None:
+        regr_set = args["regr_set"]
+    else:
+        regr_set = default_config["regr_set"]
+
+    print(regr_set)
+    config = {}
+    for rs in regr_set:
+        config[rs] = default_config["defaults"][rs]
+
+        custom_options_key = f"{rs}_options"
+        is_option = custom_options_key in args
+        if is_option and args[custom_options_key] is not None:
+            config[rs]["options"] = [args[custom_options_key]]
+        config[rs]["options"] = sort_out_options(rs, config[rs]["options"][0])
+    print(config)
+
+    dv, _, _ = wrapper(args["root_dir"], args["sub"], args["ses"], config)
+    save_dnx_vols(dv, args["odir"])

denoise/ios/vol_io.py

@@ -0,0 +1,175 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+volume manipulations
+    @author: Slava Karolis (slava.karolis@kcl.ac.uk)
+"""
+
+import numpy as np
+import nibabel as nb
+import os.path as op
+import numpy.matlib as ml
+
+
+def load_vol(path, third_var=False):
+    im_obj = nb.load(path)
+    im = im_obj.get_fdata()
+
+    dims = np.array(im.shape)
+    im = im_flatten(im)
+
+    if third_var is True:
+        return im, dims, im_obj
+    else:
+        return im, dims
+
+
+def im_flatten(im):
+    # turns volumes to time by voxels orientation, opposite to im_gather
+    if len(im.shape) == 4:
+        Nt = im.shape[3]
+        im = im.reshape((-1, Nt))
+    else:
+        im = im.reshape(-1, 1)
+    im = np.transpose(im)
+
+    return im
+
+
+def im_gather(im, dims, mask=None):
+    # turns flattened (time by voxels) data to volumetric shape
+    if len(im.shape) == 1:
+        im = im.reshape(-1, 1)
+
+    if im.shape[0] < im.shape[1]:
+        im = np.transpose(im)
+    if len(dims) == 3:
+        dims = np.append(dims, 1)
+
+    if type(mask) != type(None):
+
+        assert len(mask) == np.prod(dims[0:3]), "mask and dims are not compatable"
+        temp = np.zeros((len(mask), im.shape[1]), dtype="<f4")
+        temp[mask, :] = im
+        im = temp
+    im = im.reshape((dims))
+    return im
+
+
+def get_vol(path, vol_name=None, applymask=False):
+    # high-level addition to load_vol
+    # allow for two ways of calling volume and three ways of calling mask
+    if vol_name == None:
+        assert op.isfile(path), "Provide full path or vol_name"
+        Ac = path
+
+    else:
+        Ac = path + "/" + vol_name
+
+    path = op.split(Ac)[0]
+    #    print(Ac)
+    assert op.exists(Ac), "volume {0} doesn" "t exist".format(Ac)
+
+    func, dims = load_vol(Ac)
+
+    if isinstance(applymask, bool):
+        if applymask:
+            #  assumes that mask is in the same folder
+            assert op.exists(
+                path + "/mask.nii.gz"
+            ), "mask is not found in the same folder"
+            mask = load_vol(path + "/mask.nii.gz")[0]
+            # mask=np.array(mask)
+            mask = np.reshape(mask, -1)
+        else:
+            Nvox = func.shape[1]
+            mask = np.ones(Nvox)
+    else:
+        if type(applymask) == type("str"):  # applymask is a full path
+            assert op.exists(applymask), "mask doesn" "t exist"
+            mask = load_vol(applymask)[0]
+            mask = np.reshape(mask, -1)
+        elif type(applymask) == type(np.zeros(1)):  # applymask is  a vector
+            mask = applymask
+        else:
+            assert 1 == 0, "Check masking options"
+
+    mask = np.array(mask, dtype="bool")
+    func = func[:, mask]
+    return func, dims, mask
+
+
+def save_vol(data, path, mask=None, path_to_template=None, dims=None):
+
+    assert (type(dims) != type(None)) or (
+        type(path_to_template) != type(None)
+    ), "needs either template or dimensions of the image"
+
+    if type(path_to_template) != type(None):
+        assert op.exists(path_to_template), "Template image doesn" "t exist"
+        template = nb.load(path_to_template)
+
+        if type(dims) != type(None):
+            print("Warning: disregarding template dimensions, using dims variable")
+        else:
+            dims = template.shape
+
+    if isinstance(mask, str):
+        # mask is a full path
+        assert op.exists(mask), "mask doesn" "t exist"
+        mask = load_vol(mask)[0]
+        mask = np.reshape(mask, -1)
+    if mask is not None:
+        mask = np.array(mask, dtype="bool")
+    data = im_gather(data, dims, mask)
+
+    if path_to_template != None:
+        data = nb.Nifti1Image(data, template.affine, template.header)
+
+    else:
+        data = nb.Nifti1Image(data, np.eye(4))
+
+    nb.save(data, path)
+
+
+def select_folder(args, var_name):
+    for i in range(0, len(args)):
+        if var_name in args[i]:
+            full_path = args[i].replace(var_name, "")
+            folder = op.dirname(full_path)
+            break
+
+    return folder
+
+
+def image_name(args, var_name):
+    for i in range(0, len(args)):
+        if var_name in args[i]:
+            full_path = args[i].replace(var_name, "")
+            filez = op.basename(full_path)
+            break
+
+    basename = filez.replace(".nii.gz", "")
+    return basename
+
+
+def im_coord(dims, mask=None):
+    if len(dims) == 4:
+        dims = dims[0:3]
+
+    xyz = []
+
+    for i in range(3):
+        dims_ = np.roll(dims, -i)
+        order = np.roll([0, 1, 2], i)
+        a = np.arange(0, dims_[0]).reshape(dims_[0], 1)
+        a = ml.repmat(a, 1, dims_[1])
+        a = ml.repmat(a, 1, dims_[2])
+        x = a.reshape(dims_)
+        xf = np.transpose(x, order)
+        xyz.append(im_flatten(xf))
+
+    if mask is not None:
+        xyz = [i[:, mask] for i in xyz]
+
+    return xyz

denoise/methods/dctm.py

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Sep 11 10:33:18 2020
+
+@author: slava karolis
+"""
+
+import math
+import numpy as np
+
+
+def dctm(Ns=350, TR=2.2, lp=None, hp=None):
+    n = np.arange(0, Ns)
+
+    # check constistency
+    if lp is not None and hp is not None:
+        assert lp > hp, "Lowpass value needs to be larger than highpass"
+
+    if lp is None:
+        range_lp = []
+    else:
+        range_lp = np.arange(math.ceil(2 * (Ns * TR) / (1 / lp)), Ns)
+
+    if hp is None:
+        range_hp = []
+    else:
+        range_hp = np.arange(
+            1, math.ceil(2 * (Ns * TR) / (1 / hp))
+        )  # ceil instead of matlab's "fix"
+
+    K = np.hstack((range_hp, range_lp))
+
+    if len(K) > 0:
+        C = np.zeros((Ns, len(K)))
+        for i, k in enumerate(K):
+
+            C[:, i] = [
+                math.sqrt(2 / Ns) * math.cos(math.pi * (2 * j + 1) * (k) / (2 * Ns))
+                for j in n
+            ]
+    else:
+        C = []
+
+    return C

denoise/methods/dnx.py

@@ -0,0 +1,373 @@
+#!/usr/bin/env fslpython
+# -*- coding: utf-8 -*-
+"""
+    @author: Slava Karolis (slava.karolis@kcl.ac.uk) 
+    
+    This work is jointly copyrighted by University of Oxford and King's College London
+    Copyright 2024
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+        http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+"""
+
+import glob
+import os.path as op
+import numpy as np
+from ios.vol_io import get_vol
+import methods.mp as mp
+from methods.ic_regs import ic_regress
+
+# from scipy.io import loadmat
+from scipy.stats import zscore
+from methods.dctm import dctm
+from sklearn.decomposition import PCA
+
+
+# PATH CONSTRUCTORS
+def generate_path_to_input_folder(sub, ses, structure, generic_root):
+    # the data are expected to be organised as follows:
+    # generic_root
+    # |___ rawdata
+    # |        |___ f'sub-{sub}/ses-{ses}/func
+    # |__  derivatives
+    #          |__ dhcp_fmri_pipeline
+    #                 |__ f'sub-{sub}/ses-{ses}/func
+    assert (
+        structure == "derivatives" or structure == "rawdata" or structure is None
+    ), "structure should be rawdata or derivatives"
+    if structure == "derivatives":
+        structure = f"{structure}/dhcp_fmri_pipeline"
+    return f"{generic_root}/{structure}/sub-{sub}/ses-{ses}/func"
+
+
+def generate_path_to_dv_volume(sub, ses, generic_root):
+    input_folder = generate_path_to_input_folder(sub, ses, "rawdata", generic_root)
+    dv = glob.glob(
+        f"{input_folder}/sub-{sub}_ses-{ses}_run-*_task-rest_rec-mcdc_bold.nii.gz"
+    )
+    return dv[0]
+
+
+def generate_iv_file_name(sub, ses, name_ending, structure, generic_root=None):
+    # generates path to a regressor file
+    if structure == "rawdata":  # then need to figure out run number from the data
+        # search for the mcdc file and extract run number from that
+        dv = generate_path_to_dv_volume(sub, ses, generic_root)
+        run = op.basename(dv)
+        run = run.replace(f"sub-{sub}_ses-{ses}_run-", "")
+        run = run.replace("_task-rest_rec-mcdc_bold.nii.gz", "")
+        input_file = f"sub-{sub}_ses-{ses}_run-{run}_{name_ending}"
+    else:
+        input_file = f"sub-{sub}_ses-{ses}_{name_ending}"
+    return input_file
+
+
+def generate_path_to_mask(sub, ses, root_folder):
+    # HARDWIRED, change as appropriate.
+    # Assumes that there is a mask_dil.nii.gz file, which is a dilated version of the mask.nii.gz
+    mask_path = f"{root_folder}/derivatives/dhcp_fmri_pipeline/sub-{sub}/ses-{ses}/func/mask_dil.nii.gz"
+    assert op.exists(
+        mask_path
+    ), f"Mask file {mask_path} does not exist (not included in the release). Create it using fslmaths {root_folder}/derivatives/dhcp_fmri_pipeline/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-rest_desc-brain_mask.nii.gz -dilM {root_folder}/derivatives/dhcp_fmri_pipeline/sub-{sub}/ses-{ses}/func/mask_dil.nii.gz"
+    return mask_path
+
+
+# METHODS TO LOAD
+def load_nothing(non_arg):
+    return []
+
+
+def loc_load_cs_text(fname):
+    text_file = open(fname)
+    lines = text_file.readlines()
+    if len(lines) == 0:
+        data = None
+    else:
+        bad_vols = lines[0].split(",")
+        data = [int(i) for i in bad_vols]
+    return data
+
+
+def loc_get_vol(path, vol_name=None, applymask=False):
+    if isinstance(vol_name, list):
+        for i, name in enumerate(vol_name):
+            func_, dims, mask = get_vol(path, name, applymask)
+            if i == 0:
+                shp = list(func_.shape)
+                shp.append(len(vol_name))
+                shp = tuple(shp)
+                func = np.zeros(shp)
+            func[:, :, i] = func_
+    else:
+        func, dims, mask = get_vol(path, vol_name, applymask)
+    return func
+
+
+## IV transform methods
+def tm_bptf(iv, options, censored_vols=None):
+    # iv is a place holder
+    iv = dctm(350, 2.2, options[0], options[1])
+    return iv
+
+
+def tm_vol_censor(iv, options, censored_vols=None):
+    if options is None:
+        options = 350  # options is number of volumes here
+
+    if iv is None:
+        iv = np.zeros((options, 0))
+    else:
+        iv_ = np.zeros((options, len(iv)))
+        for i, d in enumerate(iv):
+            iv_[d, i] = 1
+        iv = iv_
+
+    return iv
+
+
+def tm_zscoring(iv, options, censored_vols=None):
+    if options is None:
+        options = iv.shape[1]
+
+    options = min(options, iv.shape[1])
+    iv = iv[:, np.arange(0, options)]
+    iv = zscore(iv, axis=0)
+    return iv
+
+
+def tm_mp_map(iv, options):
+    iv = mp.mp_expansion(iv, detrend=False, aroma=False)
+    options = min(options, iv.shape[1])
+    iv = iv[:, np.arange(0, options)]
+    iv = zscore(iv, axis=0)
+    return iv
+
+
+def do_pca(iv, thres=0.99):
+    pca = PCA(n_components=iv.shape[1] - 1, svd_solver="full")
+    pca.fit(iv)
+    a = np.cumsum(pca.explained_variance_ratio_)
+
+    if isinstance(thres, float) and thres < 1:
+        ind = [i for i, k in enumerate(a) if k > thres]
+        if len(ind) == 0:
+            Ncomps = iv.shape[1] - 1
+        else:
+            Ncomps = ind[0]
+    elif isinstance(thres, int):
+        Ncomps = thres
+    else:
+        raise "Thres needs to be either int of number of comps or < 1 (percentile)"
+
+    pc = np.transpose(pca.components_)
+    iv_demean = iv - pca.mean_
+    iv = np.dot(iv_demean, pc[:, 0:Ncomps])
+    return iv, pca
+
+
+def tm_mp_pca(iv, options, censored_vols=None):
+    Ncols = int(iv.shape[1] / 2)
+    iv1 = tm_mp_map(iv[:, 0:Ncols], options[0])
+    iv2 = tm_mp_map(iv[:, Ncols:], options[1])
+    iv = np.hstack((iv1, iv2))
+    thres = options[2]
+
+    if censored_vols is not None:
+        good_vols = np.arange(0, 350)
+        good_vols = np.setdiff1d(good_vols, censored_vols)
+
+        iv_ = np.delete(iv, censored_vols, axis=0)
+
+        iv_ = do_pca(iv_, thres)[0]
+        iv = np.zeros((iv.shape[0], iv_.shape[1]))
+        iv[good_vols, :] = iv_
+
+    else:
+        iv = do_pca(iv, thres)[0]
+
+    return iv
+
+
+def transform_function_mapping(function_name):
+    transform_func = {
+        "vol_censor": tm_vol_censor,
+        "zscoring": tm_zscoring,
+        "mp_pca": tm_mp_pca,
+        "bptf": tm_bptf,
+    }
+    return transform_func[function_name]
+
+
+####################################################################
+
+
+class IV:
+    def __init__(self, sub, ses, root_folder, config):
+        self.root_folder = root_folder
+        self.mask_path = generate_path_to_mask(sub, ses, root_folder)
+        self.sub = sub
+        self.ses = ses
+        self.regr_set = [k for k in config.keys()]
+        self.options = [config[k]["options"] for k in config.keys()]
+        self.name = [
+            config[k]["name"] for k in config.keys()
+        ]  # list of lists  for file names to load the data
+        self.structure = [
+            config[k]["structure"][0] for k in config.keys()
+        ]  # rawdata or derivatives
+        self.ismap = [config[k]["ismap"][0] for k in config.keys()]  # voxel-wise or not
+
+        for i, k in enumerate(self.options):
+            if len(k) == 1:
+                self.options[i] = k[0]
+
+        self.transform_func = []
+        for k in self.regr_set:
+            func_name = config[k]["transform_method"][0]
+            self.transform_func.append(transform_function_mapping(func_name))
+
+        self.folder_name, self.basename = self.obtain_basename()
+        self.load_func, self.load_args = self.load_method()
+
+    def obtain_basename(self):
+        folder_name = []
+        basename = []
+        for name, structure in zip(self.name, self.structure):
+            folder_name.append(
+                generate_path_to_input_folder(
+                    self.sub, self.ses, structure, self.root_folder
+                )
+            )
+            basename.append(
+                [
+                    generate_iv_file_name(
+                        self.sub, self.ses, k, structure, self.root_folder
+                    )
+                    for k in name
+                ]
+            )
+        return folder_name, basename
+
+    def load_method(self):  # defines methods to load
+        load_func = []
+        load_args = []
+        for foldname, bname, ismap, rs in zip(
+            self.folder_name, self.basename, self.ismap, self.regr_set
+        ):
+            if ismap:
+                func = loc_get_vol
+                args = {
+                    "path": foldname,
+                    "vol_name": bname,
+                    "applymask": self.mask_path,
+                }
+            else:
+                # TO CHECK HOW NONE BEHAVES FOR BFT
+                assert len(bname) == 1, "multiple or empty inputs are not allowed"
+                if rs == "vc":
+                    func = loc_load_cs_text
+                    args = {"fname": f"{foldname}/{bname[0]}"}
+                elif rs == "ica":
+                    func = np.loadtxt
+                    args = {"fname": f"{foldname}/{bname[0]}"}
+                elif rs == "bp":  # nothing to load, placeholder
+                    func = load_nothing
+                    args = {"non_arg": []}
+                else:
+                    raise NameError("no regressor is specified with this name")
+
+            load_func.append(func)
+            load_args.append(args)
+        return load_func, load_args
+
+    def loader(self):  # loads IV data
+        data = []
+        for i, func in enumerate(self.load_func):
+            print(self.load_args[i])
+            data_ = func(**self.load_args[i])
+            data.append(data_)
+        self.data = data
+
+        ind = [j for j, rs in enumerate(self.regr_set) if rs == "vc"]
+        if len(ind) == 1:
+            self.censored_vols = self.data[ind[0]]
+        else:
+            self.censored_vols = None
+
+    def transformer(self):  # ,nvox=None
+        self.transformed = [None] * len(self.data)
+
+        # populate regressors which are not voxel-specific, i.e., no need to compute them each time
+        ind = [i for i, x in enumerate(self.ismap) if x == False]
+        for i in ind:
+            self.transformed[i] = self.transform_func[i](
+                self.data[i], options=self.options[i]
+            )
+
+        if any(
+            self.ismap
+        ):  # if there are voxel specific regressors, we will will be creating regressor set dynamically
+            ind = [i for i, j in enumerate(self.ismap) if j]
+            nvox = self.data[ind[0]].shape[1]
+            self.iv = map(lambda vox: self.transformer_pervox(vox), range(nvox))
+
+        else:  # if all of them are not voxel-specific, we can freeze the regressor set
+            self.iv = np.hstack(self.transformed)
+
+    def transformer_pervox(self, vox):
+        ind = [i for i, x in enumerate(self.ismap) if x]
+        for i in ind:
+            self.transformed[i] = self.transform_func[i](
+                self.data[i][:, vox, :], self.options[i], self.censored_vols
+            )
+        return np.hstack(self.transformed)
+
+
+class Regress:
+    def __init__(self):
+        self.regress_func = ic_regress
+        self.args = {"m": None}
+
+
+class DV:
+    def __init__(self, sub, ses, generic_root):
+        self.path = generate_path_to_dv_volume(sub, ses, generic_root)
+        self.mask_path = generate_path_to_mask(sub, ses, generic_root)
+
+    def loader(self):
+        self.data, self.dims, _ = get_vol(self.path, None, self.mask_path)
+        # self.data=[self.data]
+
+    def run_regress_sequentially(self, Regress, IV):
+        resid = np.zeros(self.data.shape)
+        yhat = np.zeros(self.data.shape)
+
+        for i, iv in enumerate(IV.iv):
+            out = Regress.regress_func(iv, self.data[:, i], **Regress.args)
+            resid[:, i] = out[0].reshape(-1)
+            yhat[:, i] = out[1].reshape(-1)
+
+        return resid, yhat
+
+    def activate_dnx(self, IV, Regress):
+        self.loader()
+        IV.loader()
+        IV.transformer()
+        print("start denoise")
+        if any(IV.ismap):
+            resid, yhat = self.run_regress_sequentially(Regress, IV)
+        else:
+            out = Regress.regress_func(IV.iv, self.data, **Regress.args)
+            resid = out[0]
+            yhat = out[1]
+        self.yhat = yhat
+        self.resid = resid

denoise/methods/ic_regs.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env fslpython
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 10 19:00:08 2020
+% basic operations for temporal and spatial regressions for ICs
+@author: slava karolis
+"""
+import numpy as np
+import scipy
+from sklearn.decomposition import PCA
+from scipy.stats import t
+
+
+def regstats_part1(DM):
+    Q, R = np.linalg.qr(DM)
+    p = DM.shape[1]
+    nobs = DM.shape[0]
+    ri = np.linalg.solve(R, np.eye(p))
+    xtxi = ri @ np.transpose(ri)
+    dfe = nobs - p
+    return dfe, xtxi
+
+
+def regstats_part2(DM, beta, y, dfe, xtxi):
+    yhat = DM @ beta
+    residuals = y - yhat
+    sse = np.power(np.linalg.norm(residuals, 2), 2)
+    mse = sse / dfe
+    covb = xtxi * mse
+    se = np.sqrt(np.diag(covb))
+
+    return se
+
+
+def pinv_func(modes):
+    try:
+        pinv_DM = scipy.linalg.pinv(modes)
+    except:
+        rank = np.linalg.matrix_rank(modes[:, :-1])
+        pca = PCA(n_components=rank, svd_solver="full")
+        pca.fit(modes[:, :-1])
+        modes_tr = np.transpose(pca.components_)
+        modes = np.hstack((modes_tr, modes[:, -1].reshape(-1, 1)))
+        pinv_DM = scipy.linalg.pinv(modes)
+
+    return pinv_DM
+
+
+def ic_modes(modes, fmri):
+    func = scipy.linalg.pinv
+    if modes.shape[0] == fmri.shape[0]:
+        modes = np.c_[modes, np.ones((modes.shape[0], 1))]
+        pinv_DM = func(modes)
+        # pinv_DM=np.linalg.pinv(modes)
+        betas = pinv_DM @ fmri
+    elif modes.shape[1] == fmri.shape[1]:
+        modes = np.r_[modes, np.ones((1, modes.shape[1]))]
+        pinv_DM = func(modes)
+        # pinv_DM=np.linalg.pinv(modes)
+        betas = fmri @ pinv_DM
+    else:
+        print("matrices do not commute")
+    return betas
+
+
+def ic_regress(modes, fmri, m=None):
+    #  modes - dv to regress out from the data, may be either tc or spatial
+    #  map. Differentiated internally by their dimensionality
+    #
+    # NOTE: pinv is calculated for all modes, even if just some of modes
+    # need to be regressed out. This is to align with fsl_regfilt
+    if len(fmri.shape) == 1:
+        fmri = fmri.reshape(-1, 1)
+
+    betas = ic_modes(modes, fmri)
+
+    if m == None:
+        if modes.shape[0] == fmri.shape[0]:
+            m = np.arange(0, modes.shape[1])
+        else:
+            m = np.arange(0, modes.shape[0])
+
+    if modes.shape[0] == fmri.shape[0]:
+        # % tc are predictors, extracting spatial map
+        betas = np.delete(betas, -1, 0)
+        yhat = modes[:, m] @ betas[m, :]
+    else:
+        # % spatial maps are predictors, extracting tc
+        betas = np.delete(betas, -1, 1)
+        yhat = betas[:, m] @ modes[m, :]
+
+    fmri_resid = fmri - yhat
+
+    return fmri_resid, yhat, betas
+
+
+def regstats(IV, DV):
+
+    DM = np.c_[IV, np.ones((IV.shape[0], 1))]
+    pinv_DM = np.linalg.pinv(DM)
+    dfe, xtxi = regstats_part1(DM)
+
+    betas = pinv_DM @ DV
+    se = regstats_part2(DM, betas, DV, dfe, xtxi)
+
+    tval = [betas.tolist()[i] / se[i] for i in range(len(betas))]
+
+    pval = [t.cdf(-abs(i), dfe) * 2 for i in tval]
+
+    return tval, pval, betas

denoise/methods/mp.py

@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon May 11 17:26:18 2020
+
+@author: slava karolis
+"""
+
+
+import numpy as np
+
+
+def mp_expansion(Tn, detrend=False, aroma=False):
+    if detrend:
+        #  modelling drift as a function of time and time squared
+        DM = np.arange(Tn.shape[0]).reshape((-1, 1))
+        DM = DM - np.mean(DM, axis=0)
+        DM = np.concatenate((DM, np.power(DM, 2)), axis=1)
+        import methods.ic_regs as icr
+
+        Tn, yhat, _ = icr.ic_regress(DM, Tn)
+
+    diff_Tn = np.diff(Tn, axis=0)
+    diff_Tn = np.vstack((np.zeros((1, 6)), diff_Tn))
+
+    # center both raw and difference (to make orthogonal to squared)
+    Tn = Tn - np.mean(Tn, axis=0)
+    diff_Tn = diff_Tn - np.mean(diff_Tn, axis=0)
+
+    Tn = np.concatenate((Tn, np.power(Tn, 2), diff_Tn, np.power(diff_Tn, 2)), axis=1)
+    if aroma:
+        Tn = aroma_type(Tn)
+    else:
+        if not aroma:
+            diff_Tn = np.diff(diff_Tn, axis=0)
+            diff_Tn = np.vstack((np.zeros((1, 6)), diff_Tn))
+            diff_Tn = diff_Tn - np.mean(diff_Tn, axis=0)
+            Tn = np.concatenate((Tn, diff_Tn, np.power(diff_Tn, 2)), axis=1)
+    if detrend:
+        Tn = np.concatenate((Tn, yhat), axis=1)
+
+    return Tn
+
+
+def mp_volumetric(T, detrend=False, aroma=False):
+
+    # this guy takes per-slice mps and turnes them into per-volume via pca
+
+    # print(T.shape)
+
+    from sklearn.decomposition import PCA
+
+    pca = PCA(n_components=1)
+
+    T = np.squeeze(T)
+
+    Tn = np.zeros((T.shape[0], T.shape[2]))
+    for i in range(0, 6):
+        Tn[:, i] = pca.fit(T[:, :, i]).transform(T[:, :, i]).reshape(-1)
+
+    Tn = mp_expansion(Tn, detrend=False, aroma=aroma)
+
+    return Tn
+
+
+def aroma_type(Tn):
+
+    forward_shift = np.vstack((np.zeros((1, 24)), Tn[:-1, :24]))
+    back_shift = np.vstack((Tn[1:, :24], np.zeros((1, 24))))
+    Tn_new = np.concatenate((Tn[:, :24], forward_shift, back_shift), axis=1)
+    if Tn.shape[1] > 24:
+        Tn_new = np.concatenate((Tn_new, Tn[:, 24:]), axis=1)
+
+    return Tn_new