From 2fa93a4707d2d01c183b21331bbea246336b1ddf Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 11 Oct 2024 11:20:03 -0400 Subject: [PATCH 01/13] Initial devel notebook --- notebooks/cpca_devel.ipynb | 386 +++++++++++++++++++++++++++++++++++++ 1 file changed, 386 insertions(+) create mode 100644 notebooks/cpca_devel.ipynb diff --git a/notebooks/cpca_devel.ipynb b/notebooks/cpca_devel.ipynb new file mode 100644 index 0000000..560ea2c --- /dev/null +++ b/notebooks/cpca_devel.ipynb @@ -0,0 +1,386 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "23e7595a-8d55-4c2e-b82f-81ff4d6f158d", + "metadata": {}, + "outputs": [], + "source": [ + "input_files = '/data/SFIMJGC_Introspec/cpca_vs_introspec/spring_2024/data/cpca/JAVIER_POST_SFN2024/cpca_input_files_SUBJ_AWARE.denoise.b0.txt'\n", + "file_format = 'nifti'\n", + "n_comps = 3260\n", + "mask_fp = '/data/SFIMJGC_Introspec/pdn/PrcsData/cpca/Schaefer2018_400Parcels_7Networks_AAL2.MASK.nii.gz'\n", + "out_prefix = 'A1'\n", + "pca_type = 'complex'\n", + "rotate = None\n", + "recon = False\n", + "normalize = 'zscore'\n", + "bandpass = False\n", + "low_cut = 0.01\n", + "high_cut = 0.1\n", + "tr = None\n", + "n_bins = 30\n", + "verbose = True" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c7c67358-144e-4b0c-adc6-2ab7a8f7e7f3", + "metadata": {}, + "outputs": [], + "source": [ + "from utils.load_write import load_data, write_out\n", + "from scipy.signal import hilbert\n", + "import fbpca\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "33e577cf-4f51-4ef1-aa16-ed46abb9377a", + "metadata": {}, + "outputs": [], + "source": [ + "def hilbert_transform(input_data, verbose):\n", + " if verbose:\n", + " print('applying hilbert transform')\n", + " # hilbert transform\n", + " input_data = hilbert(input_data, axis=0)\n", + " return input_data.conj()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "699465fa-da41-4dac-b4a0-35db79d0422e", + "metadata": {}, + "outputs": [], + "source": [ + "def pca(input_data, n_comps, verbose, n_iter=10):\n", + " # compute pca\n", + " print('performing PCA/CPCA')\n", + " # get number of observations\n", + " n_samples = input_data.shape[0]\n", + " print(' number of samples = %d' % n_samples)\n", + " print(' input_data.shape[1] = %d' % input_data.shape[1])\n", + " #matrix_rank = np.linalg.matrix_rank(input_data)\n", + " #print(' rank of input matrix = % s' % str(matrix_rank))\n", + " # fbpca pca\n", + " (U, s, Va) = fbpca.pca(input_data, k=n_comps, n_iter=n_iter)\n", + " # calc explained variance\n", + " explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1]\n", + " total_var = explained_variance_.sum()\n", + " # compute PC scores\n", + " pc_scores = input_data @ Va.T\n", + " # get loadings from eigenvectors\n", + " loadings = Va.T @ np.diag(s)\n", + " loadings /= np.sqrt(input_data.shape[0]-1)\n", + " # package outputs\n", + " output_dict = {'U': U,\n", + " 's': s,\n", + " 'Va': Va,\n", + " 'loadings': loadings.T,\n", + " 'exp_var': explained_variance_,\n", + " 'pc_scores': pc_scores,\n", + " 'n_samples': n_samples,\n", + " 'n_positions': input_data.shape[1],\n", + " 'total_var': total_var}\n", + " return output_dict" + ] + }, + { + "cell_type": "markdown", + "id": "d235450c-7624-478f-8092-0ab4cffaa302", + "metadata": {}, + "source": [ + "***\n", + "### Load the Original Data" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "2cd34963-45d1-428b-b346-709c373cd719", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "initializing matrix of size (3260, 122767)\n", + "loading and concatenating 5 scans\n" + ] + } + ], + "source": [ + "func_data, mask, header = load_data(\n", + " input_files, file_format, mask_fp, normalize,\n", + " bandpass, low_cut, high_cut, tr, verbose\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "a049284b-8525-4bf2-97e3-07cf669da446", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(3260, 122767)\n" + ] + } + ], + "source": [ + "print(func_data.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "17c19249-c8f4-458b-be2d-af52acbfa4f5", + "metadata": {}, + "outputs": [], + "source": [ + "write_out(func_data,mask,header,'nifti','/data/SFIMJGC_Introspec/cpca_vs_introspec/spring_2024/data/cpca/JAVIER_POST_SFN2024/ORIG')" + ] + }, + { + "cell_type": "markdown", + "id": "d93f61df-53e7-450d-bf7f-7df89229ffa3", + "metadata": {}, + "source": [ + "***\n", + "\n", + "### Apply Hilbert Transform --> obtain analytical signal" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "09efbcef-2dc7-4f7f-95e6-177e110aae53", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "applying hilbert transform\n" + ] + } + ], + "source": [ + "# if pca_type is complex, compute hilbert transform\n", + "if pca_type == 'complex':\n", + " func_data = hilbert_transform(func_data, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "e5294ef1-baf1-42ae-8eb1-33076199cacd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(3260, 122767)" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "func_data.shape" + ] + }, + { + "cell_type": "markdown", + "id": "c5c143f8-6b1c-40d3-b78d-397a39b39ac0", + "metadata": {}, + "source": [ + "***\n", + "\n", + "### Apply PCA" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "e1ba7b70-9b22-4e48-8505-ab06cbb60e21", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "performing PCA/CPCA\n", + " number of samples = 3260\n", + " input_data.shape[1] = 122767\n" + ] + } + ], + "source": [ + "pca_output = pca(func_data, n_comps, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "338607a1-4cc8-49a2-a705-4d7245202ae8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((3260, 3260), (3260,), (3260, 122767))" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pca_output['U'].shape, pca_output['s'].shape, pca_output['Va'].shape" + ] + }, + { + "cell_type": "markdown", + "id": "9e120bcc-e14c-45e9-8f27-cbbbfe575c0b", + "metadata": {}, + "source": [ + "***\n", + "### Remove Components" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "a0c4b0a3-f25a-45fd-a1fb-1adc98722ef4", + "metadata": {}, + "outputs": [], + "source": [ + "pca_output['s_mod'] = pca_output['s'].copy()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "53ca3ef0-6cdf-4022-a685-1b176058f7f8", + "metadata": {}, + "outputs": [], + "source": [ + "pca_output['s_mod'][:3] = 0" + ] + }, + { + "cell_type": "markdown", + "id": "3712a3c7-3bb3-4f4d-966b-9ecd96123fba", + "metadata": {}, + "source": [ + "***\n", + "### Reconstruct the data" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "f209406b-c28d-470e-a772-9801798061d4", + "metadata": {}, + "outputs": [], + "source": [ + "func_data_mod_reconstructed = np.dot(pca_output['U'] * pca_output['s_mod'], pca_output['Va'])" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "469b0f5b-ed43-4b77-abe3-b2d906a092a4", + "metadata": {}, + "outputs": [], + "source": [ + "func_data_modified = np.real(func_data_mod_reconstructed)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "62afbfaa-f733-44a2-abbe-da53af9ee1c9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[-4.08930025e-14, 2.49407012e-14, 3.17258923e-14, ...,\n", + " 0.00000000e+00, 4.27470732e-02, 0.00000000e+00],\n", + " [ 1.06886709e-13, 4.12176498e-14, 4.84417944e-14, ...,\n", + " 0.00000000e+00, -2.46599703e-02, 0.00000000e+00],\n", + " [-3.61867827e-14, -6.14628260e-14, -8.21240862e-15, ...,\n", + " 0.00000000e+00, -7.14048055e-02, 0.00000000e+00],\n", + " ...,\n", + " [ 4.11187276e-14, -1.50444896e-13, 1.22669481e-14, ...,\n", + " 0.00000000e+00, -9.09199494e-01, 0.00000000e+00],\n", + " [ 1.45857796e-14, -9.76858223e-14, -2.39312872e-14, ...,\n", + " 0.00000000e+00, -7.04478752e-01, 0.00000000e+00],\n", + " [ 7.01746098e-14, -4.18684331e-14, 2.98126921e-14, ...,\n", + " 0.00000000e+00, -1.93273426e-01, 0.00000000e+00]])" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "func_data_modified" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "07edeca7-e0e7-4974-b7d6-24422ddef6b2", + "metadata": {}, + "outputs": [], + "source": [ + "write_out(func_data_modified,mask,header,'nifti','/data/SFIMJGC_Introspec/cpca_vs_introspec/spring_2024/data/cpca/JAVIER_POST_SFN2024/RECON_3OUT.nii')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31d4e7d8-8b10-408b-9e87-d86c71c8e4d0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "BOLD WAVES 2024a", + "language": "python", + "name": "bold_waves_2024a" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From fbafa8808e8f4e9c717a32495349377278820387 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 11 Oct 2024 11:50:56 -0400 Subject: [PATCH 02/13] read_input_file now supports several entries per line --- utils/load_write.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/utils/load_write.py b/utils/load_write.py index 60abac0..50a7bbb 100644 --- a/utils/load_write.py +++ b/utils/load_write.py @@ -189,9 +189,11 @@ def read_input_file(input_files): fps = [line.rstrip() for line in file] # remove extra lines, if any fps = [line for line in fps if len(line)>0] + # Only separate into input output if there is more than one entry per line + if '\t' in fps[0]: + fps = [item.split('\t') for item in fps] return fps - def write_out(data, mask, header, file_format, out_prefix): # write out brain maps to nifti or cifti if file_format == 'nifti': From 8a2943eb95f4b3bf2b93d592e9e6e29d0c9fed2c Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 11 Oct 2024 13:12:18 -0400 Subject: [PATCH 03/13] load_data now can handle having only input paths or both input and ouput paths up to load_scans --- utils/load_write.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/utils/load_write.py b/utils/load_write.py index 50a7bbb..7bbf2f0 100644 --- a/utils/load_write.py +++ b/utils/load_write.py @@ -61,13 +61,23 @@ def load_data(input_files, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose): # read file paths from input .txt file fps = read_input_file(input_files) + # master function for loading and concatenating functional scans # parameters check parameter_check(fps, file_format, tr, bandpass, mask_fp) + # Separate input and ouput paths + if isinstance(fps[0],list): + fps_inputs = [i[0] for i in fps] + fps_outputs = [i[1] for i in fps] + else: + fps_inputs = fps + fps_outputs = None + del fps + # Pull file paths data, mask, header = load_scans( - fps, file_format, mask_fp, normalize, bandpass, + fps_inputs, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose ) From a3a17a0e7dc7eede32293a43cfbfcdac51fa3531 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 11 Oct 2024 13:23:32 -0400 Subject: [PATCH 04/13] data_trs now available --- utils/load_write.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/utils/load_write.py b/utils/load_write.py index 7bbf2f0..b1a66e2 100644 --- a/utils/load_write.py +++ b/utils/load_write.py @@ -76,16 +76,18 @@ def load_data(input_files, file_format, mask_fp, normalize, del fps # Pull file paths - data, mask, header = load_scans( + data, mask, header, data_trs = load_scans( fps_inputs, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose ) - return data, mask, header + return data, mask, header, data_trs def load_scans(fps, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose): + # scan duration in number of trs + scan_trs = [] # get # of scans n_scans = len(fps) # if file_format = 'nifti', load mask @@ -115,10 +117,13 @@ def load_scans(fps, file_format, mask_fp, normalize, tr, verbose) # get # observations data_n = data.shape[0] + scan_trs.append(data_n) # Normalize data before concatenation if normalize == 'zscore': + print(f' zscoring {fp}') data = zscore(data, nan_policy='omit') elif normalize == 'mean_center': + print(f' mean centering {fp}') data = data - np.mean(data, axis=0) # fill nans w/ 0 in regions of poor functional scan coverage data = np.nan_to_num(data) @@ -127,7 +132,7 @@ def load_scans(fps, file_format, mask_fp, normalize, # increase counter by # of observations in subject scan indx += data_n - return group_data, mask, header + return group_data, mask, header, scan_trs def load_file(fp, file_format, mask, bandpass, From 32321c4e55d7c94a241c6deb8b3774ce600eb386 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Wed, 23 Oct 2024 09:02:48 -0400 Subject: [PATCH 05/13] functions now accept multiple input paths (in and out) --- utils/load_write.py | 79 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 65 insertions(+), 14 deletions(-) diff --git a/utils/load_write.py b/utils/load_write.py index b1a66e2..fa024e9 100644 --- a/utils/load_write.py +++ b/utils/load_write.py @@ -7,7 +7,7 @@ from scipy.io import loadmat from scipy.stats import zscore from scipy.signal import butter, sosfiltfilt - +import os.path as osp def butter_bandpass(lowcut, highcut, fs, order=5): # create butterworth bandpass filter based on lowcut and highcut @@ -59,29 +59,48 @@ def initialize_matrix(fps, nscans, file_format, mask, verbose): def load_data(input_files, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose): + """This function loads resting-state data and concatenates it in time. + + Inputs + ====== + + Outputs + ======= + data: concatenated data + mask: mask + header: header + data_trs: number of acquisitions per scan + input_paths: path to input files + out_asis_paths: path to output files (all components) + out_removed_paths: path to output files (after removal of components) + """ # read file paths from input .txt file fps = read_input_file(input_files) - - # master function for loading and concatenating functional scans - # parameters check - parameter_check(fps, file_format, tr, bandpass, mask_fp) - + input_paths, out_asis_paths, out_removed_paths = None, None, None # Separate input and ouput paths if isinstance(fps[0],list): - fps_inputs = [i[0] for i in fps] - fps_outputs = [i[1] for i in fps] + print(' + [load_data]: both input and output paths were provided') + if len(fps[0]) != 3: + print(' ++ ERROR: Not enought entries per line in input file.') + exit(1) + input_paths = [i[0] for i in fps] + out_asis_paths = [i[1] for i in fps] + out_removed_paths = [i[2] for i in fps] else: - fps_inputs = fps - fps_outputs = None + input_paths = fps del fps + # master function for loading and concatenating functional scans + # parameters check + parameter_check(input_paths, file_format, tr, bandpass, mask_fp) + # Pull file paths data, mask, header, data_trs = load_scans( - fps_inputs, file_format, mask_fp, normalize, bandpass, + input_paths, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose ) - return data, mask, header, data_trs + return data, mask, header, data_trs, input_paths, out_asis_paths, out_removed_paths def load_scans(fps, file_format, mask_fp, normalize, @@ -205,8 +224,13 @@ def read_input_file(input_files): # remove extra lines, if any fps = [line for line in fps if len(line)>0] # Only separate into input output if there is more than one entry per line + print(' + [read_input_file]: Number of entries in input file: %d lines' % (len(fps))) if '\t' in fps[0]: fps = [item.split('\t') for item in fps] + print(' + [read_input_file]: Number of paths per input: %d paths' % (len(fps[0]))) + elif ' ' in fps[0]: + fps = [item.split(' ') for item in fps] + print(' + [read_input_file]: Number of paths per input: %d paths' % (len(fps[0]))) return fps def write_out(data, mask, header, file_format, out_prefix): @@ -219,7 +243,9 @@ def write_out(data, mask, header, file_format, out_prefix): nifti_4d[mask_bin, :] = data.T # write out brain w/ nibabel and mask affine nifti_out = nb.Nifti2Image(nifti_4d, mask.affine) - nb.save(nifti_out, f'{out_prefix}.nii') + if (not out_prefix.endswith('.nii')) & (not out_prefix.endswith('.nii.gz')): + out_prefix = out_prefix + '.nii' + nb.save(nifti_out, out_prefix) elif file_format == 'cifti': # https://neurostars.org/t/alter-size-of-matrix-for-new-cifti-header-nibabel/20903/2 # first get tr from Series axis @@ -233,4 +259,29 @@ def write_out(data, mask, header, file_format, out_prefix): # need to create a new header due to change in matrix shape nb.save(cifti_out, f'{out_prefix}.dtseries.nii') - +def write_modified_scans(data,mask,header,file_format,data_trs,file_paths, verbose): + ## # extracting output paths + ## fps = read_input_file(file_paths) + ## # Separate input and ouput paths + ## if isinstance(fps[0],list): + ## fps_inputs = [i[0] for i in fps] + ## fps_outputs = [i[1] for i in fps] + ## else: + ## fps_inputs = fps + ## fps_outputs = None + ## del fps + + # initlializing tr index + indx = 0 + n_files_to_write = len(file_paths) + for i,out_path in enumerate(file_paths): + ii = i + 1 + out_data = data[indx:(indx+data_trs[i]),:] + if verbose: + print(f'[{ii}/{n_files_to_write}] writing modified data {out_path}') + out_dir = osp.dirname(out_path) + if not osp.exists(out_dir): + print(f'[{ii}/{n_files_to_write}] --> Creating output folder first: {out_dir}') + os.makedirs(out_dir) + write_out(out_data,mask,header,file_format,out_path) + indx += data_trs[i] From b63e26ec70fba7c919053c7b18f5458c5c42ef3e Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Wed, 23 Oct 2024 09:04:02 -0400 Subject: [PATCH 06/13] new functions to save decomposition results as hdf5 object (too big to save as matlab) --- cpca.py | 160 +++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 136 insertions(+), 24 deletions(-) diff --git a/cpca.py b/cpca.py index cb55946..7e9c206 100644 --- a/cpca.py +++ b/cpca.py @@ -1,20 +1,40 @@ import argparse import numpy as np +import pandas as pd import fbpca import warnings - +import h5py from numpy.linalg import pinv from scipy.io import savemat from scipy.signal import hilbert from scipy.stats import zscore from utils.cpca_reconstruction import cpca_recon -from utils.load_write import load_data, write_out +from utils.load_write import load_data, write_out, write_modified_scans from xmca.tools.rotation import varimax, promax +# def save_dict_to_hdf5(filename, data_dict): +# with h5py.File(filename, 'w') as h5file: +# for key, value in data_dict.items(): +# h5file.create_dataset(key, data=value) + +def save_dict_to_hdf5(filename, data_dict): + with h5py.File(filename, 'w') as h5file: + for key, value in data_dict.items(): + if isinstance(value, np.ndarray): + # Check if the array contains objects + if value.dtype == np.dtype('O'): + # Try to convert the object array to a numeric array if possible + try: + value = np.array(value, dtype=np.float64) + except ValueError: + print(f"Skipping key {key} - cannot convert to numeric.") + continue + h5file.create_dataset(key, data=value) + else: + print(f"Skipping key {key} - not a numpy array.") + def hilbert_transform(input_data, verbose): - if verbose: - print('applying hilbert transform') # hilbert transform input_data = hilbert(input_data, axis=0) return input_data.conj() @@ -49,31 +69,46 @@ def package_parameters(n_comps, mask_fp, file_format, return params -def pca(input_data, n_comps, verbose, n_iter=10): +def pca(input_data, n_comps, pca_type, verbose, n_iter=10): # compute pca print('performing PCA/CPCA') # get number of observations - n_samples = input_data.shape[0] - print(' number of samples = %d' % n_samples) - print(' input_data.shape[1] = %d' % input_data.shape[1]) + n_samples = input_data.shape[0] + n_vertices = input_data.shape[1] + print(' number of samples = %d' % n_samples) + print(' number of vertices/voxels = %d' % n_vertices) #matrix_rank = np.linalg.matrix_rank(input_data) #print(' rank of input matrix = % s' % str(matrix_rank)) # fbpca pca (U, s, Va) = fbpca.pca(input_data, k=n_comps, n_iter=n_iter) # calc explained variance - explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1] + # 1. Get eigs + eigs = (s ** 2) / (n_samples-1) + # 2. Compute Variance Explained from eigenvaluesa + # Assumptions: + # a) We are working with complex data --> we count the number of vertices twice + # b) Input were normalized timseries --> each vertex/voxel has a variance = 1 + if pca_type == 'complex': + explained_variance_ = np.array([eig/(n_vertices*2) for eig in eigs]) + if pca_type == 'real': + explained_variance_ = np.array([eig/(n_vertices) for eig in eigs]) + # Original Formulation + # explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1] total_var = explained_variance_.sum() - # compute PC scores + + # 3. Compute PC scores pc_scores = input_data @ Va.T # get loadings from eigenvectors loadings = Va.T @ np.diag(s) loadings /= np.sqrt(input_data.shape[0]-1) - # package outputs + + # 4. Package outputs output_dict = {'U': U, 's': s, 'Va': Va, 'loadings': loadings.T, 'exp_var': explained_variance_, + 'eigs': eigs, 'pc_scores': pc_scores, 'n_samples': n_samples, 'n_positions': input_data.shape[1], @@ -100,7 +135,7 @@ def rotation(pca_output, data, rotation, verbose): def write_results(pca_output, pca_type, mask, file_format, - header, rotate, out_prefix): + header, rotate, out_prefix, n_comps_to_save): # write out results of pca analysis # create output name if out_prefix is None if out_prefix is None: @@ -112,8 +147,23 @@ def write_results(pca_output, pca_type, mask, file_format, out_prefix += f'_{rotate}' out_prefix += '_results' + # Write variance explained + out_df = None + if 's_modified' in pca_output: + out_df = pd.DataFrame(np.vstack([pca_output['s'], pca_output['eigs'], pca_output['exp_var'], pca_output['s_modified']]).T,columns=['s','eigs','exp_var','s_modified']) + else: + out_df = pd.DataFrame(np.vstack([pca_output['s'], pca_output['eigs'], pca_output['exp_var']]).T,columns=['s','eigs','exp_var']) + out_df.index.name = 'cpca_id' + out_df_path = f'{out_prefix}_exp_var.txt' + out_df.to_csv(out_df_path) + print(" + [write_results]: Explained Variance Dataframe written to disk [%s]" % out_df_path) + # get loadings - loadings = pca_output['loadings'] + loadings = pca_output['loadings'][0:n_comps_to_save:] + print('=====================') + print(loadings.shape) + print(type(pca_output)) + print('=====================') # Write brain maps if file_format in ('nifti', 'cifti'): if pca_type == 'complex': @@ -139,17 +189,25 @@ def write_results(pca_output, pca_type, mask, file_format, loadings, mask, header, file_format, out_prefix ) # write out pca results dictionary to .mat file - savemat(f'{out_prefix}.mat', pca_output) - + print(" + Starting to save pca_output as h5 object",end=' ') + save_dict_to_hdf5(f'{out_prefix}.h5', pca_output) + #savemat(f'{out_prefix}.mat', pca_output) + print("[DONE]") def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, pca_type, rotate, recon, normalize, bandpass, - low_cut, high_cut, tr, n_bins, verbose): + low_cut, high_cut, tr, n_bins, verbose,recon_data,n_comps_to_remove, n_comps_to_recon,save_pca_out): print('++ Entering Run cpca...') + print(' + number of components to recon = %s' % str(n_comps_to_recon)) + print(' + number of components to compute = %s' % str(n_comps)) + print(' + number of components to remove = %s' % str(n_comps_to_remove)) + print(' + recon data after component removal? %s' % str(recon_data)) + print(' + recon components separately? %s' % str(recon)) print(' + bandpass = %s' % str(bandpass)) print(' + verbose = %s' % str(verbose)) # load dataset - func_data, mask, header = load_data( + print("++ Loading data into memory.....") + func_data, mask, header, func_data_trs, input_paths, out_asis_paths, out_removed_paths = load_data( input_files, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose ) @@ -158,19 +216,50 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, print(' + Applying Hilbert Transform ...') func_data = hilbert_transform(func_data, verbose) + # if n_comps not provided, set it to maximum possible + if n_comps is None: + n_comps = np.min(func_data.shape) + print(f' + Automatically setting n_comps = {n_comps}') # compute pca - pca_output = pca(func_data, n_comps, verbose) + pca_output = pca(func_data, n_comps, pca_type, verbose) # rotate pca weights, if specified if rotate is not None: print(' + Applying rotation ...') pca_output = rotation(pca_output, func_data, rotate, verbose) + # free memory + print(' + Freeing memory by deleting the func_data variable') + del func_data + + # if cpca, and recon_data = True, create reconstructed data after removal of complex PC components + # ================================================================================================ + if recon_data & (pca_type == 'complex') & (n_comps_to_remove is not None): + if verbose: + print(f'performing data reconstruction after removal of {n_comps_to_remove} components') + # Nulling the first n_comps_to_remove components + pca_output['s_modified'] = pca_output['s'].copy() + pca_output['s_modified'][:n_comps_to_remove] = 0 + # Reconstructing the data (in analytical form) + func_data_modified = np.dot(pca_output['U'] * pca_output['s_modified'], pca_output['Va']) + func_data_allcomps = np.dot(pca_output['U'] * pca_output['s'], pca_output['Va']) + # Extracting the real part of the reconstructed data + func_data_modified = np.real(func_data_modified) + func_data_allcomps = np.real(func_data_allcomps) + # Write reconstructed data to disk + write_modified_scans(func_data_modified,mask,header,file_format,func_data_trs, out_removed_paths, verbose) + write_modified_scans(func_data_allcomps,mask,header,file_format,func_data_trs, out_asis_paths, verbose) + # free memory + print(' + Freeing memory by deleting the func_data_modified variable') + del func_data_modified + del func_data_allcomps + # if cpca, and recon=True, create reconstructed time courses of complex PC + # ======================================================================== if recon & (pca_type == 'complex'): if verbose: print('performing CPCA component time series reconstruction') - if n_comps > 10: + if n_comps_to_recon > 10: warnings.warn( 'the # of components estimated is large, CPCA reconstruction ' 'will create a separate file for each component. This may take ' @@ -178,7 +267,7 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, ) del func_data # free up memory cpca_recon(pca_output, rotate, file_format, - mask, header, out_prefix, n_bins) + mask, header, out_prefix, n_bins, n_comps_to_recon) elif recon & (pca_type == 'real'): warnings.warn('Time series reconstruction only available for CPCA') @@ -190,8 +279,10 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, # write out results if verbose: print('writing out results') - write_results(pca_output, pca_type, mask, file_format, - header, rotate, out_prefix) + if save_pca_out: + print(' writing pca_output dictorionary to disk') + write_results(pca_output, pca_type, mask, file_format, + header, rotate, out_prefix, n_comps_to_recon) if __name__ == '__main__': @@ -206,7 +297,18 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, type=str) parser.add_argument('-n', '--n_comps', help=' Number of components from PCA', - required=True, + required=False, + default=None, + type=int) + parser.add_argument('-n_comps_to_recon', '--n_comps_to_recon', + help=' Number of components from PCA to reconstruct (one at a time)', + required=False, + default=None, + type=int) + parser.add_argument('-rn', '--n_comps_to_remove', + help=' Number of components to remove when reconstructing the data', + required=False, + default=None, type=int) parser.add_argument('-m', '--mask', help='path to brain mask in nifti format. Only needed ' @@ -243,6 +345,14 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, 'complex PCA', action='store_true', required=False) + parser.add_argument('-recon_data', '--recon_data', + help='Whether to reconstruct the data after removal of compoments', + action='store_true', + required=False) + parser.add_argument('-save_pca_out', '--save_pca_out', + help='Save the PCA dictionary', + action='store_true', + required=False) parser.add_argument('-norm', '--normalize', help='Type of scan normalization before group ' 'concatenation. It is recommend to z-score', @@ -285,6 +395,7 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, help='turn off printing', action='store_false', required=False) + args_dict = vars(parser.parse_args()) run_cpca(args_dict['input'], args_dict['n_comps'], args_dict['mask'], @@ -293,4 +404,5 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, args_dict['recon'], args_dict['normalize'], args_dict['bandpass_filter'], args_dict['bandpass_filter_low'], args_dict['bandpass_filter_high'], args_dict['sampling_unit'], - args_dict['n_recon_bins'], args_dict['verbose_off']) + args_dict['n_recon_bins'], args_dict['verbose_off'],args_dict['recon_data'],args_dict['n_comps_to_remove'], + args_dict['n_comps_to_recon'],args_dict['save_pca_out']) From 744338c6d17a230b3b712cae907cc025f04bfa11 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Wed, 23 Oct 2024 09:04:37 -0400 Subject: [PATCH 07/13] cpca_recon now accept option to choose how many components (less than the total) to recon --- utils/cpca_reconstruction.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/cpca_reconstruction.py b/utils/cpca_reconstruction.py index db7d05e..abbc851 100644 --- a/utils/cpca_reconstruction.py +++ b/utils/cpca_reconstruction.py @@ -58,11 +58,11 @@ def write_results(phase_ts, rotation, n, mask, header, def cpca_recon(cpca_res, rotation, file_format, mask, header, - out_prefix, n_bins): + out_prefix, n_bins, n_comps_to_recon): # reconstruct cpca component 'movies' from cpca results bin_indx_all = [] bin_centers_all = [] - n_comps = cpca_res['pc_scores'].shape[1] + n_comps = n_comps_to_recon #cpca_res['pc_scores'].shape[1] for n in range(n_comps): recon_ts = reconstruct_ts(cpca_res, [n], rotation) phase_ts = np.angle(cpca_res['pc_scores'][:,n]) From bd73f0165b97dd69a2a5ce2e04dde770696be1b6 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Wed, 23 Oct 2024 10:46:45 -0400 Subject: [PATCH 08/13] load_write while zscoring now says on which scan out of how many it is working --- utils/load_write.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/load_write.py b/utils/load_write.py index fa024e9..069739e 100644 --- a/utils/load_write.py +++ b/utils/load_write.py @@ -129,7 +129,7 @@ def load_scans(fps, file_format, mask_fp, normalize, # initialize counter indx=0 # Loop through files and concatenate/append - for fp in fps: + for n,fp in enumerate(fps): # load file data, header = load_file(fp, file_format, mask_bin, bandpass, low_cut, high_cut, @@ -139,7 +139,7 @@ def load_scans(fps, file_format, mask_fp, normalize, scan_trs.append(data_n) # Normalize data before concatenation if normalize == 'zscore': - print(f' zscoring {fp}') + print(f' [{n}/{n_scans}] zscoring {fp}') data = zscore(data, nan_policy='omit') elif normalize == 'mean_center': print(f' mean centering {fp}') From 619add508c510e76604ee9fcb8ce3b9e6ccfd21b Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Wed, 23 Oct 2024 10:47:06 -0400 Subject: [PATCH 09/13] var explained now reported as percentage --- cpca.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/cpca.py b/cpca.py index 7e9c206..a49ea14 100644 --- a/cpca.py +++ b/cpca.py @@ -89,9 +89,9 @@ def pca(input_data, n_comps, pca_type, verbose, n_iter=10): # a) We are working with complex data --> we count the number of vertices twice # b) Input were normalized timseries --> each vertex/voxel has a variance = 1 if pca_type == 'complex': - explained_variance_ = np.array([eig/(n_vertices*2) for eig in eigs]) + explained_variance_ = 100*np.array([eig/(n_vertices*2) for eig in eigs]) if pca_type == 'real': - explained_variance_ = np.array([eig/(n_vertices) for eig in eigs]) + explained_variance_ = 100*np.array([eig/(n_vertices) for eig in eigs]) # Original Formulation # explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1] total_var = explained_variance_.sum() @@ -196,7 +196,7 @@ def write_results(pca_output, pca_type, mask, file_format, def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, pca_type, rotate, recon, normalize, bandpass, - low_cut, high_cut, tr, n_bins, verbose,recon_data,n_comps_to_remove, n_comps_to_recon,save_pca_out): + low_cut, high_cut, tr, n_bins, verbose,recon_data,n_comps_to_remove, n_comps_to_recon,save_pca_out, calc_rank): print('++ Entering Run cpca...') print(' + number of components to recon = %s' % str(n_comps_to_recon)) print(' + number of components to compute = %s' % str(n_comps)) @@ -205,6 +205,7 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, print(' + recon components separately? %s' % str(recon)) print(' + bandpass = %s' % str(bandpass)) print(' + verbose = %s' % str(verbose)) + print(' + calculate rank of data = %s' % str(calc_rank)) # load dataset print("++ Loading data into memory.....") func_data, mask, header, func_data_trs, input_paths, out_asis_paths, out_removed_paths = load_data( @@ -215,11 +216,20 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, if pca_type == 'complex': print(' + Applying Hilbert Transform ...') func_data = hilbert_transform(func_data, verbose) - + # if requested, calculate the rank of the data + if calc_rank == True: + func_data_rank = np.linalg.matrix_rank(func_data) + print('f + Input data rank = %d' % func_data_rank) + # if n_comps not provided, set it to maximum possible if n_comps is None: - n_comps = np.min(func_data.shape) - print(f' + Automatically setting n_comps = {n_comps}') + if calc_rank == True: + n_comps = func_data_rank + print(f' + Automatically setting n_comps = rank ==> n_comps = {n_comps}') + else: + n_comps = np.min(func_data.shape) + print(f' + Automatically setting n_comps = min(func_data.shape) ==> n_comps = {n_comps}') + # compute pca pca_output = pca(func_data, n_comps, pca_type, verbose) @@ -365,6 +375,10 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, 'a butterworth filter', action='store_true', required=False) + parser.add_argument('-calc_rank', '--calc_rank', + help='Whether to estimate the rank of the input data', + action='store_true', + required=False) parser.add_argument('-f_low', '--bandpass_filter_low', help='Low cut frequency for bandpass filter in Hz', required=False, @@ -405,4 +419,4 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, args_dict['bandpass_filter'], args_dict['bandpass_filter_low'], args_dict['bandpass_filter_high'], args_dict['sampling_unit'], args_dict['n_recon_bins'], args_dict['verbose_off'],args_dict['recon_data'],args_dict['n_comps_to_remove'], - args_dict['n_comps_to_recon'],args_dict['save_pca_out']) + args_dict['n_comps_to_recon'],args_dict['save_pca_out'],args_dict['calc_rank']) From b47eeb90f5ecd588b5d26c82e13cf6ae7491816a Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 17 Jan 2025 12:48:50 -0500 Subject: [PATCH 10/13] cpca now computes explained variance stlightly differently --- cpca.py | 99 +- notebooks/cpca_devel.ipynb | 1995 +++++++++++++++++++++++++++++++++++- utils/load_write.py | 13 +- 3 files changed, 2015 insertions(+), 92 deletions(-) diff --git a/cpca.py b/cpca.py index a49ea14..662280a 100644 --- a/cpca.py +++ b/cpca.py @@ -4,15 +4,15 @@ import fbpca import warnings import h5py +import os #CW added from numpy.linalg import pinv from scipy.io import savemat from scipy.signal import hilbert -from scipy.stats import zscore +from scipy.stats import zscore, tvar from utils.cpca_reconstruction import cpca_recon from utils.load_write import load_data, write_out, write_modified_scans from xmca.tools.rotation import varimax, promax - # def save_dict_to_hdf5(filename, data_dict): # with h5py.File(filename, 'w') as h5file: # for key, value in data_dict.items(): @@ -69,32 +69,35 @@ def package_parameters(n_comps, mask_fp, file_format, return params -def pca(input_data, n_comps, pca_type, verbose, n_iter=10): +def pca(input_data,total_var_in_func_data, n_comps, pca_type, verbose, n_iter=2): # compute pca - print('performing PCA/CPCA') + print(' + [pca]: Entereing the PCA function --> performing PCA/CPCA') # get number of observations n_samples = input_data.shape[0] n_vertices = input_data.shape[1] - print(' number of samples = %d' % n_samples) - print(' number of vertices/voxels = %d' % n_vertices) - #matrix_rank = np.linalg.matrix_rank(input_data) - #print(' rank of input matrix = % s' % str(matrix_rank)) + print(' + [pca] INFO: number of samples = %d' % n_samples) + print(' + [pca] INFO: number of vertices/voxels = %d' % n_vertices) + print(' + [pca] INFO: Total variance in the data = %.f' % total_var_in_func_data) # fbpca pca (U, s, Va) = fbpca.pca(input_data, k=n_comps, n_iter=n_iter) # calc explained variance # 1. Get eigs eigs = (s ** 2) / (n_samples-1) - # 2. Compute Variance Explained from eigenvaluesa + # 2. Compute Variance Explained from eigenvalues + explained_variance_ = 100*np.array([eig/total_var_in_func_data for eig in eigs]) + total_var = explained_variance_.sum() + print(' + [pca] INFO: Total variance explained in PCA = %s %%' % str(total_var)) + + # OLDER CODE THAT IS INNACURATE WHEN YOU HAVE VOXELS WITH FLAT TIMESERIES # Assumptions: # a) We are working with complex data --> we count the number of vertices twice # b) Input were normalized timseries --> each vertex/voxel has a variance = 1 - if pca_type == 'complex': - explained_variance_ = 100*np.array([eig/(n_vertices*2) for eig in eigs]) - if pca_type == 'real': - explained_variance_ = 100*np.array([eig/(n_vertices) for eig in eigs]) + #if pca_type == 'complex': + # explained_variance_ = 100*np.array([eig/(n_vertices*2) for eig in eigs]) + #if pca_type == 'real': + # explained_variance_ = 100*np.array([eig/(n_vertices) for eig in eigs]) # Original Formulation # explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1] - total_var = explained_variance_.sum() # 3. Compute PC scores pc_scores = input_data @ Va.T @@ -197,56 +200,78 @@ def write_results(pca_output, pca_type, mask, file_format, def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, pca_type, rotate, recon, normalize, bandpass, low_cut, high_cut, tr, n_bins, verbose,recon_data,n_comps_to_remove, n_comps_to_recon,save_pca_out, calc_rank): - print('++ Entering Run cpca...') - print(' + number of components to recon = %s' % str(n_comps_to_recon)) - print(' + number of components to compute = %s' % str(n_comps)) - print(' + number of components to remove = %s' % str(n_comps_to_remove)) - print(' + recon data after component removal? %s' % str(recon_data)) - print(' + recon components separately? %s' % str(recon)) - print(' + bandpass = %s' % str(bandpass)) - print(' + verbose = %s' % str(verbose)) - print(' + calculate rank of data = %s' % str(calc_rank)) + print('++ [run_cpca]: Entering Run cpca...') + print(' + number of components to recon = %s' % str(n_comps_to_recon)) + print(' + number of components to compute = %s' % str(n_comps)) + print(' + number of components to remove = %s' % str(n_comps_to_remove)) + print(' + recon data after component removal? %s' % str(recon_data)) + print(' + recon components separately? %s' % str(recon)) + print(' + bandpass = %s' % str(bandpass)) + print(' + verbose = %s' % str(verbose)) + print(' + calculate rank of data = %s' % str(calc_rank)) # load dataset - print("++ Loading data into memory.....") + print(" + [run_cpca]: Loading data into memory.....") func_data, mask, header, func_data_trs, input_paths, out_asis_paths, out_removed_paths = load_data( input_files, file_format, mask_fp, normalize, bandpass, low_cut, high_cut, tr, verbose ) + # if pca_type is complex, compute hilbert transform if pca_type == 'complex': - print(' + Applying Hilbert Transform ...') + print(' + [run_cpca]: Applying Hilbert Transform ...') func_data = hilbert_transform(func_data, verbose) + + # CW: save hilbert transformed data to confirm same with new environment + #func_dict = {} + #func_dict['func_data'] = func_data + #save_dict_to_hdf5(f'{out_prefix}_func_data.h5',func_dict) + + # Estimate and write total variance in the input data (once concatenated and hilbert) + #var_in_func_data = func_data.var(axis=0) + print(' + [run_cpca]: finished Hilbert Transform') + #var_in_func_data = tvar(func_data, axis=0, ddof=0) # CW added + var_in_func_data = np.nanvar(func_data, axis=0) # CW added + print(' + [run_cpca]: calculated variance') + total_var_in_func_data = var_in_func_data.sum() + print(' + [run_cpca]: summed variance') + total_var_in_func_data_df = pd.DataFrame(var_in_func_data,columns=['Voxelwise Variance']) + total_var_in_func_data_df.index.name = 'voxel_ID' + + total_var_in_func_data_df_path = f'{out_prefix}_{pca_type}_total_var_in_func_data.txt' + total_var_in_func_data_df.to_csv(total_var_in_func_data_df_path) + print(" + [run_cpca]: Wrote variance of original data in [%s]" % total_var_in_func_data_df_path) + # if requested, calculate the rank of the data if calc_rank == True: func_data_rank = np.linalg.matrix_rank(func_data) - print('f + Input data rank = %d' % func_data_rank) + print(' + [run_cpca]: Input data rank = %d' % func_data_rank) # if n_comps not provided, set it to maximum possible if n_comps is None: if calc_rank == True: n_comps = func_data_rank - print(f' + Automatically setting n_comps = rank ==> n_comps = {n_comps}') + print(f' + [run_cpca]: Automatically setting n_comps = rank ==> n_comps = {n_comps}') else: n_comps = np.min(func_data.shape) - print(f' + Automatically setting n_comps = min(func_data.shape) ==> n_comps = {n_comps}') + print(f' + [run_cpca]: Automatically setting n_comps = min(func_data.shape) ==> n_comps = {n_comps}') # compute pca - pca_output = pca(func_data, n_comps, pca_type, verbose) + pca_output = pca(func_data, total_var_in_func_data, n_comps, pca_type, verbose, n_iter=8) # rotate pca weights, if specified if rotate is not None: - print(' + Applying rotation ...') + print(' + [run_cpca]: Applying rotation ...') pca_output = rotation(pca_output, func_data, rotate, verbose) # free memory - print(' + Freeing memory by deleting the func_data variable') + print(' + [run_cpca]: Freeing memory by deleting the func_data variable') del func_data # if cpca, and recon_data = True, create reconstructed data after removal of complex PC components # ================================================================================================ if recon_data & (pca_type == 'complex') & (n_comps_to_remove is not None): if verbose: - print(f'performing data reconstruction after removal of {n_comps_to_remove} components') + print(f' + [run_cpca]: Performing data reconstruction after removal of {n_comps_to_remove} components') # Nulling the first n_comps_to_remove components pca_output['s_modified'] = pca_output['s'].copy() pca_output['s_modified'][:n_comps_to_remove] = 0 @@ -260,7 +285,7 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, write_modified_scans(func_data_modified,mask,header,file_format,func_data_trs, out_removed_paths, verbose) write_modified_scans(func_data_allcomps,mask,header,file_format,func_data_trs, out_asis_paths, verbose) # free memory - print(' + Freeing memory by deleting the func_data_modified variable') + print(' + [run_cpca]: Freeing memory by deleting the func_data_modified variable') del func_data_modified del func_data_allcomps @@ -268,7 +293,7 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, # ======================================================================== if recon & (pca_type == 'complex'): if verbose: - print('performing CPCA component time series reconstruction') + print(' + [run_cpca]: performing CPCA component time series reconstruction') if n_comps_to_recon > 10: warnings.warn( 'the # of components estimated is large, CPCA reconstruction ' @@ -288,13 +313,12 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, ) # write out results if verbose: - print('writing out results') + print(' + [run_cpca]: writing out results') if save_pca_out: - print(' writing pca_output dictorionary to disk') + print(' + writing pca_output dictorionary to disk') write_results(pca_output, pca_type, mask, file_format, header, rotate, out_prefix, n_comps_to_recon) - if __name__ == '__main__': """Run complex or standard principal component analysis""" parser = argparse.ArgumentParser(description='Run CPCA or PCA analysis') @@ -412,6 +436,7 @@ def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, args_dict = vars(parser.parse_args()) + os.environ["MKL_INTERFACE_LAYER"] = "ILP64" #CW added run_cpca(args_dict['input'], args_dict['n_comps'], args_dict['mask'], args_dict['file_format'], args_dict['output_prefix'], args_dict['pca_type'], args_dict['rotate'], diff --git a/notebooks/cpca_devel.ipynb b/notebooks/cpca_devel.ipynb index 560ea2c..35268e5 100644 --- a/notebooks/cpca_devel.ipynb +++ b/notebooks/cpca_devel.ipynb @@ -3,13 +3,42 @@ { "cell_type": "code", "execution_count": 1, + "id": "8fb7427f-7279-455c-96fe-660ee6b25969", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.append('../')" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c7c67358-144e-4b0c-adc6-2ab7a8f7e7f3", + "metadata": {}, + "outputs": [], + "source": [ + "from utils.load_write import load_data, write_out, read_input_file\n", + "from cpca import hilbert_transform, pca\n", + "#from scipy.signal import hilbert\n", + "import fbpca\n", + "import numpy as np\n", + "import nibabel as nb\n", + "from scipy.stats import zscore\n", + "import os\n", + "import os.path as osp" + ] + }, + { + "cell_type": "code", + "execution_count": 3, "id": "23e7595a-8d55-4c2e-b82f-81ff4d6f158d", "metadata": {}, "outputs": [], "source": [ - "input_files = '/data/SFIMJGC_Introspec/cpca_vs_introspec/spring_2024/data/cpca/JAVIER_POST_SFN2024/cpca_input_files_SUBJ_AWARE.denoise.b0.txt'\n", + "input_files = '/data/SFIMJGC_Introspec/cpca_vs_introspec/post_sfn_2024/config/cpca_input_output_15.txt'\n", "file_format = 'nifti'\n", - "n_comps = 3260\n", + "n_comps = None\n", "mask_fp = '/data/SFIMJGC_Introspec/pdn/PrcsData/cpca/Schaefer2018_400Parcels_7Networks_AAL2.MASK.nii.gz'\n", "out_prefix = 'A1'\n", "pca_type = 'complex'\n", @@ -21,74 +50,939 @@ "high_cut = 0.1\n", "tr = None\n", "n_bins = 30\n", - "verbose = True" + "verbose = True\n", + "n_comps_to_remove = 3\n", + "recon_data = True" ] }, { "cell_type": "code", - "execution_count": 2, - "id": "c7c67358-144e-4b0c-adc6-2ab7a8f7e7f3", + "execution_count": 4, + "id": "f53e06f7-d236-43ca-974d-f193199376b4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "++ Entering Run cpca...\n", + " + bandpass = False\n", + " + verbose = True\n" + ] + } + ], + "source": [ + "print('++ Entering Run cpca...')\n", + "print(' + bandpass = %s' % str(bandpass))\n", + "print(' + verbose = %s' % str(verbose)) " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6c2a8020-c8c3-4fb9-b25c-67852f723773", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "initializing matrix of size (9780, 122767)\n", + "loading and concatenating 15 scans\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010024/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010072/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010005/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010159/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010231/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010233/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010103/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010115/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010060/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010137/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010124/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010200/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010056/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010067/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010142/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n" + ] + } + ], + "source": [ + "# load dataset\n", + "func_data, mask, header, func_data_trs = load_data(\n", + " input_files, file_format, mask_fp, normalize, \n", + " bandpass, low_cut, high_cut, tr, verbose\n", + ") " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "94adc160-9b9e-40bc-a1be-779c75dc5204", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " + Applying Hilbert Transform ...\n", + "CPU times: user 1min 57s, sys: 15.3 s, total: 2min 12s\n", + "Wall time: 2min 12s\n" + ] + } + ], + "source": [ + "%%time\n", + "# if pca_type is complex, compute hilbert transform\n", + "if pca_type == 'complex':\n", + " print(' + Applying Hilbert Transform ...')\n", + " func_data = hilbert_transform(func_data, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "1f900a23-ead9-454b-a172-ea798be4f9a8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " + Automatically setting n_comps = 9780\n", + "performing PCA/CPCA\n", + " number of samples = 9780\n", + " input_data.shape[1] = 122767\n" + ] + }, + { + "ename": "ValueError", + "evalue": "failed in converting hidden `rwork' of _flapack.zgesdd to C/Fortran array", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "\u001b[0;31mValueError\u001b[0m: failed to create intent(cache|hide)|optional array -- must have defined dimensions, but dims[0] = -1702338196", + "\nThe above exception was the direct cause of the following exception:\n", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "File \u001b[0;32m:6\u001b[0m\n", + "File \u001b[0;32m/gpfs/gsfs11/users/SFIMJGC_Introspec/cpca_vs_introspec/sw/complex_pca/notebooks/../cpca.py:60\u001b[0m, in \u001b[0;36mpca\u001b[0;34m(input_data, n_comps, verbose, n_iter)\u001b[0m\n\u001b[1;32m 56\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m input_data.shape[1] = \u001b[39m\u001b[38;5;132;01m%d\u001b[39;00m\u001b[38;5;124m'\u001b[39m \u001b[38;5;241m%\u001b[39m input_data\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m1\u001b[39m])\n\u001b[1;32m 57\u001b[0m \u001b[38;5;66;03m#matrix_rank = np.linalg.matrix_rank(input_data)\u001b[39;00m\n\u001b[1;32m 58\u001b[0m \u001b[38;5;66;03m#print(' rank of input matrix = % s' % str(matrix_rank))\u001b[39;00m\n\u001b[1;32m 59\u001b[0m \u001b[38;5;66;03m# fbpca pca\u001b[39;00m\n\u001b[0;32m---> 60\u001b[0m (U, s, Va) \u001b[38;5;241m=\u001b[39m \u001b[43mfbpca\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpca\u001b[49m\u001b[43m(\u001b[49m\u001b[43minput_data\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mk\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mn_comps\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mn_iter\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mn_iter\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 61\u001b[0m \u001b[38;5;66;03m# calc explained variance\u001b[39;00m\n\u001b[1;32m 62\u001b[0m explained_variance_ \u001b[38;5;241m=\u001b[39m ((s \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m \u001b[38;5;241m2\u001b[39m) \u001b[38;5;241m/\u001b[39m (n_samples \u001b[38;5;241m-\u001b[39m \u001b[38;5;241m1\u001b[39m)) \u001b[38;5;241m/\u001b[39m input_data\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m1\u001b[39m]\n", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/fbpca.py:1644\u001b[0m, in \u001b[0;36mpca\u001b[0;34m(A, k, raw, n_iter, l)\u001b[0m\n\u001b[1;32m 1640\u001b[0m \u001b[38;5;66;03m#\u001b[39;00m\n\u001b[1;32m 1641\u001b[0m \u001b[38;5;66;03m# SVD the centered A directly if l >= m/1.25 or l >= n/1.25.\u001b[39;00m\n\u001b[1;32m 1642\u001b[0m \u001b[38;5;66;03m#\u001b[39;00m\n\u001b[1;32m 1643\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m l \u001b[38;5;241m>\u001b[39m\u001b[38;5;241m=\u001b[39m m \u001b[38;5;241m/\u001b[39m \u001b[38;5;241m1.25\u001b[39m \u001b[38;5;129;01mor\u001b[39;00m l \u001b[38;5;241m>\u001b[39m\u001b[38;5;241m=\u001b[39m n \u001b[38;5;241m/\u001b[39m \u001b[38;5;241m1.25\u001b[39m:\n\u001b[0;32m-> 1644\u001b[0m (U, s, Va) \u001b[38;5;241m=\u001b[39m \u001b[43msvd\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[43mA\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtodense\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mif\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43missparse\u001b[49m\u001b[43m(\u001b[49m\u001b[43mA\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1645\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01melse\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mA\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mones\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[43mm\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mc\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfull_matrices\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m 1646\u001b[0m \u001b[38;5;66;03m#\u001b[39;00m\n\u001b[1;32m 1647\u001b[0m \u001b[38;5;66;03m# Retain only the leftmost k columns of U, the uppermost\u001b[39;00m\n\u001b[1;32m 1648\u001b[0m \u001b[38;5;66;03m# k rows of Va, and the first k entries of s.\u001b[39;00m\n\u001b[1;32m 1649\u001b[0m \u001b[38;5;66;03m#\u001b[39;00m\n\u001b[1;32m 1650\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m U[:, :k], s[:k], Va[:k, :]\n", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/scipy/linalg/_decomp_svd.py:141\u001b[0m, in \u001b[0;36msvd\u001b[0;34m(a, full_matrices, compute_uv, overwrite_a, check_finite, lapack_driver)\u001b[0m\n\u001b[1;32m 137\u001b[0m lwork \u001b[38;5;241m=\u001b[39m _compute_lwork(gesXd_lwork, a1\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m0\u001b[39m], a1\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m1\u001b[39m],\n\u001b[1;32m 138\u001b[0m compute_uv\u001b[38;5;241m=\u001b[39mcompute_uv, full_matrices\u001b[38;5;241m=\u001b[39mfull_matrices)\n\u001b[1;32m 140\u001b[0m \u001b[38;5;66;03m# perform decomposition\u001b[39;00m\n\u001b[0;32m--> 141\u001b[0m u, s, v, info \u001b[38;5;241m=\u001b[39m \u001b[43mgesXd\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma1\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcompute_uv\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcompute_uv\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlwork\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mlwork\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 142\u001b[0m \u001b[43m \u001b[49m\u001b[43mfull_matrices\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfull_matrices\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43moverwrite_a\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moverwrite_a\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 144\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m info \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[1;32m 145\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m LinAlgError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSVD did not converge\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", + "\u001b[0;31mValueError\u001b[0m: failed in converting hidden `rwork' of _flapack.zgesdd to C/Fortran array" + ] + } + ], + "source": [ + "%%time\n", + "# if n_comps not provided, set it to maximum possible\n", + "if n_comps is None:\n", + " n_comps = np.min(func_data.shape)\n", + " print(f' + Automatically setting n_comps = {n_comps}')\n", + "# compute pca\n", + "pca_output = pca(func_data, n_comps, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "47eafb2d-6db5-4b05-a6ec-db00abe5bab5", "metadata": {}, "outputs": [], "source": [ - "from utils.load_write import load_data, write_out\n", - "from scipy.signal import hilbert\n", - "import fbpca\n", - "import numpy as np" + "from scipy.linalg import svd" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 9, + "id": "b6e47dda-9f1b-4a91-a365-9b4b21a6f3c5", + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "failed in converting hidden `rwork' of _flapack.zgesdd to C/Fortran array", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "\u001b[0;31mValueError\u001b[0m: failed to create intent(cache|hide)|optional array -- must have defined dimensions, but dims[0] = -1702338196", + "\nThe above exception was the direct cause of the following exception:\n", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[9], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m U, s, Va \u001b[38;5;241m=\u001b[39m \u001b[43msvd\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc_data\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfull_matrices\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/scipy/linalg/_decomp_svd.py:141\u001b[0m, in \u001b[0;36msvd\u001b[0;34m(a, full_matrices, compute_uv, overwrite_a, check_finite, lapack_driver)\u001b[0m\n\u001b[1;32m 137\u001b[0m lwork \u001b[38;5;241m=\u001b[39m _compute_lwork(gesXd_lwork, a1\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m0\u001b[39m], a1\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m1\u001b[39m],\n\u001b[1;32m 138\u001b[0m compute_uv\u001b[38;5;241m=\u001b[39mcompute_uv, full_matrices\u001b[38;5;241m=\u001b[39mfull_matrices)\n\u001b[1;32m 140\u001b[0m \u001b[38;5;66;03m# perform decomposition\u001b[39;00m\n\u001b[0;32m--> 141\u001b[0m u, s, v, info \u001b[38;5;241m=\u001b[39m \u001b[43mgesXd\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma1\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcompute_uv\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcompute_uv\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlwork\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mlwork\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 142\u001b[0m \u001b[43m \u001b[49m\u001b[43mfull_matrices\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfull_matrices\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43moverwrite_a\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43moverwrite_a\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 144\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m info \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[1;32m 145\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m LinAlgError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSVD did not converge\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", + "\u001b[0;31mValueError\u001b[0m: failed in converting hidden `rwork' of _flapack.zgesdd to C/Fortran array" + ] + } + ], + "source": [ + "U, s, Va = svd(func_data, full_matrices=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48353aa9-966c-45cd-ac3f-37a84fe21afd", + "metadata": {}, + "outputs": [], + "source": [ + "import tensorflow as tf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16144be4-cb8b-4cf7-85ca-f5a59e220573", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cefa71e-9f4f-4ce4-baa6-99baf492f866", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c1b8d90-51db-4bbd-b0cd-b213ac356373", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c452b59-1e67-40b2-a90f-52c08a533cb2", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "4210b779-429b-4ce9-b49f-36c657bb46d8", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0e98f04-3a3b-4385-a860-b679441f67e6", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "d303b858-53c4-4087-8a4b-e84b235d6e31", + "metadata": {}, + "outputs": [], + "source": [ + "from sklearn.decomposition import IncrementalPCA, PCA" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "29acc0cd-f38f-465d-88be-dedaa29f5dbd", + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "Complex data not supported\n[[ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n 4.70789052e-16+0.1852642j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n 7.41056841e-17+0.09085743j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n 1.04619789e-16+0.05789111j 0.00000000e+00-0.j ]\n ...\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n -8.97182539e-01-0.34183792j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n -6.71154534e-01+0.2499411j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n -1.98113753e-01+0.37533229j 0.00000000e+00-0.j ]]\n", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mComplexWarning\u001b[0m Traceback (most recent call last)", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/sklearn/utils/validation.py:856\u001b[0m, in \u001b[0;36mcheck_array\u001b[0;34m(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, estimator, input_name)\u001b[0m\n\u001b[1;32m 855\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 856\u001b[0m array \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43marray\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43morder\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43morder\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 857\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m ComplexWarning \u001b[38;5;28;01mas\u001b[39;00m complex_warning:\n", + "\u001b[0;31mComplexWarning\u001b[0m: Casting complex values to real discards the imaginary part", + "\nThe above exception was the direct cause of the following exception:\n", + "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[14], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m pca \u001b[38;5;241m=\u001b[39m PCA(n_components\u001b[38;5;241m=\u001b[39mn_comps, svd_solver\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mrandomized\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m----> 2\u001b[0m X_pca \u001b[38;5;241m=\u001b[39m \u001b[43mpca\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfit\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc_data\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/sklearn/decomposition/_pca.py:408\u001b[0m, in \u001b[0;36mPCA.fit\u001b[0;34m(self, X, y)\u001b[0m\n\u001b[1;32m 385\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Fit the model with X.\u001b[39;00m\n\u001b[1;32m 386\u001b[0m \n\u001b[1;32m 387\u001b[0m \u001b[38;5;124;03mParameters\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 399\u001b[0m \u001b[38;5;124;03m Returns the instance itself.\u001b[39;00m\n\u001b[1;32m 400\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 401\u001b[0m check_scalar(\n\u001b[1;32m 402\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn_oversamples,\n\u001b[1;32m 403\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mn_oversamples\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m 404\u001b[0m min_val\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m,\n\u001b[1;32m 405\u001b[0m target_type\u001b[38;5;241m=\u001b[39mnumbers\u001b[38;5;241m.\u001b[39mIntegral,\n\u001b[1;32m 406\u001b[0m )\n\u001b[0;32m--> 408\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_fit\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 409\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\n", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/sklearn/decomposition/_pca.py:456\u001b[0m, in \u001b[0;36mPCA._fit\u001b[0;34m(self, X)\u001b[0m\n\u001b[1;32m 450\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m issparse(X):\n\u001b[1;32m 451\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(\n\u001b[1;32m 452\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPCA does not support sparse input. See \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 453\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mTruncatedSVD for a possible alternative.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 454\u001b[0m )\n\u001b[0;32m--> 456\u001b[0m X \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_validate_data\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 457\u001b[0m \u001b[43m \u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfloat64\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfloat32\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mensure_2d\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcopy\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcopy\u001b[49m\n\u001b[1;32m 458\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 460\u001b[0m \u001b[38;5;66;03m# Handle n_components==None\u001b[39;00m\n\u001b[1;32m 461\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn_components \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/sklearn/base.py:577\u001b[0m, in \u001b[0;36mBaseEstimator._validate_data\u001b[0;34m(self, X, y, reset, validate_separately, **check_params)\u001b[0m\n\u001b[1;32m 575\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mValidation should be done on X, y or both.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 576\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m no_val_X \u001b[38;5;129;01mand\u001b[39;00m no_val_y:\n\u001b[0;32m--> 577\u001b[0m X \u001b[38;5;241m=\u001b[39m \u001b[43mcheck_array\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minput_name\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mX\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mcheck_params\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 578\u001b[0m out \u001b[38;5;241m=\u001b[39m X\n\u001b[1;32m 579\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m no_val_X \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m no_val_y:\n", + "File \u001b[0;32m/data/SFIMJGC_HCP7T/Apps/envs/bold_waves_2024a/lib/python3.9/site-packages/sklearn/utils/validation.py:858\u001b[0m, in \u001b[0;36mcheck_array\u001b[0;34m(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, estimator, input_name)\u001b[0m\n\u001b[1;32m 856\u001b[0m array \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39masarray(array, order\u001b[38;5;241m=\u001b[39morder, dtype\u001b[38;5;241m=\u001b[39mdtype)\n\u001b[1;32m 857\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m ComplexWarning \u001b[38;5;28;01mas\u001b[39;00m complex_warning:\n\u001b[0;32m--> 858\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m 859\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mComplex data not supported\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(array)\n\u001b[1;32m 860\u001b[0m ) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mcomplex_warning\u001b[39;00m\n\u001b[1;32m 862\u001b[0m \u001b[38;5;66;03m# It is possible that the np.array(..) gave no warning. This happens\u001b[39;00m\n\u001b[1;32m 863\u001b[0m \u001b[38;5;66;03m# when no dtype conversion happened, for example dtype = None. The\u001b[39;00m\n\u001b[1;32m 864\u001b[0m \u001b[38;5;66;03m# result is that np.array(..) produces an array of complex dtype\u001b[39;00m\n\u001b[1;32m 865\u001b[0m \u001b[38;5;66;03m# and we need to catch and raise exception for such cases.\u001b[39;00m\n\u001b[1;32m 866\u001b[0m _ensure_no_complex_data(array)\n", + "\u001b[0;31mValueError\u001b[0m: Complex data not supported\n[[ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n 4.70789052e-16+0.1852642j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n 7.41056841e-17+0.09085743j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n 1.04619789e-16+0.05789111j 0.00000000e+00-0.j ]\n ...\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n -8.97182539e-01-0.34183792j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n -6.71154534e-01+0.2499411j 0.00000000e+00-0.j ]\n [ 0.00000000e+00-0.j 0.00000000e+00-0.j\n 0.00000000e+00-0.j ... 0.00000000e+00-0.j\n -1.98113753e-01+0.37533229j 0.00000000e+00-0.j ]]\n" + ] + } + ], + "source": [ + "pca = PCA(n_components=n_comps, svd_solver='randomized')\n", + "X_pca = pca.fit(func_data)" + ] + }, + { + "cell_type": "markdown", + "id": "9383ec30-c8dc-4221-a3fb-31bdc027ae90", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "fa8ed668-0974-4d0e-a1ca-84456e007b12", + "metadata": {}, + "outputs": [], + "source": [ + "# rotate pca weights, if specified\n", + "if rotate is not None:\n", + " print(' + Applying rotation ...')\n", + " pca_output = rotation(pca_output, func_data, rotate, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "39713e8c-7f6b-4114-9297-2d176dbda643", + "metadata": {}, + "outputs": [], + "source": [ + "def write_modified_scans(data,mask,header,file_format,data_trs,file_paths, verbose):\n", + " # extracting output paths\n", + " fps = read_input_file(file_paths)\n", + " # Separate input and ouput paths\n", + " if isinstance(fps[0],list):\n", + " fps_inputs = [i[0] for i in fps]\n", + " fps_outputs = [i[1] for i in fps]\n", + " else:\n", + " fps_inputs = fps\n", + " fps_outputs = None\n", + " del fps\n", + " \n", + " # initlializing tr index\n", + " indx = 0\n", + " n_files_to_write = len(fps_outputs)\n", + " for i,out_path in enumerate(fps_outputs):\n", + " ii = i + 1\n", + " out_data = data[indx:(indx+data_trs[i]),:]\n", + " if verbose:\n", + " print(f'[{ii}/{n_files_to_write}] writing modified data {out_path}')\n", + " out_dir = osp.dirname(out_path)\n", + " if not osp.exists(out_dir):\n", + " print(f'[{ii}/{n_files_to_write}] --> Creating output folder first: {out_dir}')\n", + " os.makedirs(out_dir)\n", + " write_out(out_data,mask,header,file_format,out_path) \n", + " indx += data_trs[i]" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "0f878277-109c-4492-a78c-286833574921", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "performing data reconstruction after removal of 3 components\n", + "[1/5] writing modified data /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010024/preprocessed/func/pb07_cpca/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.3cPCA_rem.nii.gz\n", + "[2/5] writing modified data /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010072/preprocessed/func/pb07_cpca/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.3cPCA_rem.nii.gz\n", + "[3/5] writing modified data /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010005/preprocessed/func/pb07_cpca/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.3cPCA_rem.nii.gz\n", + "[4/5] writing modified data /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010159/preprocessed/func/pb07_cpca/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.3cPCA_rem.nii.gz\n", + "[5/5] writing modified data /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010231/preprocessed/func/pb07_cpca/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.3cPCA_rem.nii.gz\n" + ] + } + ], + "source": [ + "# if cpca, and recon_data = True, create reconstructed data after removal of complex PC components\n", + "if recon_data & (pca_type == 'complex') & (n_comps_to_remove > 0):\n", + " if verbose:\n", + " print(f'performing data reconstruction after removal of {n_comps_to_remove} components')\n", + " # Nulling the first n_comps_to_remove components\n", + " pca_output['s_modified'] = pca_output['s'].copy()\n", + " pca_output['s_modified'][:n_comps_to_remove] = 0\n", + " # Reconstructing the data (in analytical form)\n", + " func_data_modified = np.dot(pca_output['U'] * pca_output['s_modified'], pca_output['Va'])\n", + " # Extracting the real part of the reconstructed data\n", + " func_data_modified = np.real(func_data_modified)\n", + " # Write reconstructed data to disk\n", + " write_modified_scans(func_data_modified,mask,header,file_format,func_data_trs, input_files, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "838eee3a-826f-490e-8c71-74ae3ef7514d", + "metadata": {}, + "outputs": [], + "source": [ + "# if cpca, and recon=True, create reconstructed time courses of complex PC\n", + "if recon & (pca_type == 'complex'):\n", + " if verbose:\n", + " print('performing CPCA component time series reconstruction')\n", + " if n_comps > 10:\n", + " warnings.warn(\n", + " 'the # of components estimated is large, CPCA reconstruction '\n", + " 'will create a separate file for each component. This may take '\n", + " 'a while.'\n", + " )\n", + " del func_data # free up memory\n", + " cpca_recon(pca_output, rotate, file_format,\n", + " mask, header, out_prefix, n_bins)\n", + "elif recon & (pca_type == 'real'):\n", + " warnings.warn('Time series reconstruction only available for CPCA')\n", + "\n", + "# put input parameters into PCA results dicitonary\n", + "pca_output['params'] = package_parameters(\n", + " n_comps, mask_fp, file_format, pca_type, rotate, \n", + " normalize, bandpass, low_cut, high_cut, tr\n", + ")\n", + "# write out results\n", + "if verbose:\n", + " print('writing out results')\n", + "write_results(pca_output, pca_type, mask, file_format, \n", + " header, rotate, out_prefix)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "42cbbf76-e5f0-43e6-a8fa-8c54ba1ab6f5", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a990c80c-cb67-4a17-a474-491671577930", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c215c06-17ed-4d0e-b433-21b9091bbcb5", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1150f277-2eb9-4eaa-b12b-61a29879eadf", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "475837cc-f398-4512-b2e9-9f40302c60a3", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ffb4d0a-bbe6-4a0d-af58-ca5ce251a433", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ab321f6-ed10-46d0-9f72-34154eaaaa67", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b786698-5e1e-4300-a183-58264059890e", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3e1ea4e-eadc-49b6-a2eb-c705526889e6", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19cd26e6-85de-47d2-838a-9bfb4b6af290", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35604e7f-e861-49f7-bbe8-524014952dd5", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5faf7679-e44f-4d03-ab5a-26cd607e9654", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "02b68829-c595-49bf-9e31-57ecc046f9e1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "initializing matrix of size (3260, 122767)\n", + "loading and concatenating 5 scans\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010024/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010072/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010005/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010159/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010231/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n" + ] + } + ], + "source": [ + "# load dataset\n", + "func_data, mask, header, func_data_trs = load_data(\n", + " input_files, file_format, mask_fp, normalize, \n", + " bandpass, low_cut, high_cut, tr, verbose) " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "b58a8ade-41cc-4249-abf5-eb60f591aabe", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "+ automatically setting n_comps = 3260\n" + ] + } + ], + "source": [ + "# if n_comps not provided, set it to maximum possible\n", + "if n_comps is None:\n", + " n_comps = np.min(func_data.shape)\n", + " print(f'+ automatically setting n_comps = {n_comps}')" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "5bae67fc-e0b0-4679-af79-c6de2a903103", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " + Applying Hilbert Transform ...\n" + ] + } + ], + "source": [ + "# if pca_type is complex, compute hilbert transform\n", + "if pca_type == 'complex':\n", + " print(' + Applying Hilbert Transform ...')\n", + " func_data = hilbert_transform(func_data, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "898e9bdf-322b-4979-b481-ff629fb5bd33", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "performing PCA/CPCA\n", + " number of samples = 3260\n", + " input_data.shape[1] = 122767\n" + ] + } + ], + "source": [ + "# compute pca\n", + "pca_output = pca(func_data, n_comps, verbose)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d5da180-8eaa-459b-8d19-46e7a1842813", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57cb04bc-4001-4acf-a373-b2ff54974e1e", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0540f59-8006-402d-a837-2d00ad6cc840", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a06b3bf2-34f1-4d53-a00c-5b129c08bd4e", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 91, "id": "33e577cf-4f51-4ef1-aa16-ed46abb9377a", "metadata": {}, "outputs": [], "source": [ - "def hilbert_transform(input_data, verbose):\n", - " if verbose:\n", - " print('applying hilbert transform')\n", - " # hilbert transform\n", - " input_data = hilbert(input_data, axis=0)\n", - " return input_data.conj()" + "def hilbert_transform(input_data, verbose):\n", + " if verbose:\n", + " print('applying hilbert transform')\n", + " # hilbert transform\n", + " input_data = hilbert(input_data, axis=0)\n", + " return input_data.conj()" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "id": "699465fa-da41-4dac-b4a0-35db79d0422e", + "metadata": {}, + "outputs": [], + "source": [ + "def pca(input_data, n_comps, verbose, n_iter=10):\n", + " # compute pca\n", + " print('performing PCA/CPCA')\n", + " # get number of observations\n", + " n_samples = input_data.shape[0]\n", + " print(' number of samples = %d' % n_samples)\n", + " print(' input_data.shape[1] = %d' % input_data.shape[1])\n", + " #matrix_rank = np.linalg.matrix_rank(input_data)\n", + " #print(' rank of input matrix = % s' % str(matrix_rank))\n", + " # fbpca pca\n", + " (U, s, Va) = fbpca.pca(input_data, k=n_comps, n_iter=n_iter)\n", + " # calc explained variance\n", + " explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1]\n", + " total_var = explained_variance_.sum()\n", + " # compute PC scores\n", + " pc_scores = input_data @ Va.T\n", + " # get loadings from eigenvectors\n", + " loadings = Va.T @ np.diag(s)\n", + " loadings /= np.sqrt(input_data.shape[0]-1)\n", + " # package outputs\n", + " output_dict = {'U': U,\n", + " 's': s,\n", + " 'Va': Va,\n", + " 'loadings': loadings.T,\n", + " 'exp_var': explained_variance_,\n", + " 'pc_scores': pc_scores,\n", + " 'n_samples': n_samples,\n", + " 'n_positions': input_data.shape[1],\n", + " 'total_var': total_var}\n", + " return output_dict" + ] + }, + { + "cell_type": "code", + "execution_count": 93, + "id": "11394257-92ad-4924-81ac-ca57dd4d5855", + "metadata": {}, + "outputs": [], + "source": [ + "def read_input_file(input_files):\n", + " # load input .txt file specifying file paths\n", + " with open(input_files, 'r') as file:\n", + " fps = [line.rstrip() for line in file]\n", + " # remove extra lines, if any\n", + " fps = [line for line in fps if len(line)>0]\n", + " # Only separate into input output if there is more than one entry per line\n", + " if '\\t' in fps[0]:\n", + " fps = [item.split('\\t') for item in fps]\n", + " return fps" + ] + }, + { + "cell_type": "code", + "execution_count": 94, + "id": "ad473a5f-c971-456a-affc-b3cd61534c45", + "metadata": {}, + "outputs": [], + "source": [ + "def parameter_check(fps, file_format, tr, bandpass, mask):\n", + " # check errors that may have not been caught by argparse\n", + " # ensure a mask is supplied if file_format='nifi'\n", + " if (mask is None) and (file_format == 'nifti'):\n", + " raise Exception(\n", + " 'a mask .nii file must be supplied when file format is nifti'\n", + " )\n", + " # ensure that TR is provided if bandpass = True\n", + " if bandpass and (tr is None):\n", + " raise Exception(\n", + " 'the TR must be supplied if bandpass=True'\n", + " )\n", + " # test whether the files in the input match those specified in file format\n", + " try:\n", + " data_obj = nb.load(fps[0])\n", + " # check whether it is nifti\n", + " if isinstance(data_obj, nb.nifti1.Nifti1Image):\n", + " actual_format = 'nifti'\n", + " elif isinstance(data_obj, nb.cifti2.cifti2.Cifti2Image):\n", + " actual_format = 'cifti'\n", + " except ImageFileError:\n", + " actual_format = 'txt'\n", + "\n", + " if file_format != actual_format:\n", + " raise Exception(\n", + " f\"\"\"\n", + " It looks like the file format specified: '{file_format}''\n", + " does not match the file format of the input files: '{actual_format}'\n", + " \"\"\"\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 96, + "id": "6079a931-857d-4ba0-a4d8-2117a9dca162", + "metadata": {}, + "outputs": [], + "source": [ + "def initialize_matrix(fps, nscans, file_format, mask, verbose):\n", + " # initialize empty matrix for faster concatenation of scans\n", + " # if file format is nifti, get # of voxels in mask\n", + " if file_format == 'nifti':\n", + " n_ts = len(np.nonzero(mask)[0])\n", + " # if file format is cifti, get # of vertices and voxels from first scan\n", + " elif file_format == 'cifti':\n", + " cifti = nb.load(fps[0])\n", + " n_ts = cifti.shape[1]\n", + " # if file format is text, get # of time courses from first scan\n", + " elif file_format == 'txt':\n", + " txt = np.loadtxt(fps[0])\n", + " n_ts = txt.shape[1]\n", + " # loop through file paths, load header and get # of observations\n", + " # and add up across files to get total # observations\n", + " n_t = 0\n", + " for fp in fps:\n", + " if file_format == 'nifti':\n", + " nifti = nb.load(fp)\n", + " n_t += nifti.header['dim'][4]\n", + " elif file_format == 'cifti':\n", + " cifti = nb.load(fp)\n", + " n_t += cifti.shape[0]\n", + " elif file_format == 'txt':\n", + " n_t += sum(1 for _ in open(fp))\n", + " if verbose:\n", + " print(f'initializing matrix of size ({n_t}, {n_ts})')\n", + " # initialize matrix with zeros\n", + " matrix_init = np.zeros((n_t, n_ts))\n", + " return matrix_init" + ] + }, + { + "cell_type": "code", + "execution_count": 98, + "id": "51628106-266e-490b-9911-cc7236111709", + "metadata": {}, + "outputs": [], + "source": [ + "def load_file(fp, file_format, mask, bandpass, \n", + " low_cut, high_cut, tr, verbose):\n", + " # Load file based on file format\n", + " if file_format == 'nifti':\n", + " nifti = nb.load(fp)\n", + " nifti_data = nifti.get_fdata()\n", + " data = convert_2d_nifti(mask, nifti_data)\n", + " header = nifti.header\n", + " elif file_format == 'cifti':\n", + " cifti = nb.load(fp)\n", + " data = cifti.get_fdata()\n", + " header = cifti.header\n", + " elif file_format == 'txt':\n", + " data = np.loadtxt(fp)\n", + " header = None \n", + "\n", + " # if bandpass = True, bandpass all time courses by (low_cut - high_cut Hz)\n", + " if bandpass:\n", + " npad=1000 # pad 1000 samples\n", + " fs = 1/tr #tr to sampling rate\n", + " sos = butter_bandpass(low_cut, high_cut, fs)\n", + " # Median padding to reduce edge effects\n", + " data_pad = np.pad(data,[(npad, npad), (0, 0)], 'median')\n", + " # backward and forward filtering\n", + " data_filt = sosfiltfilt(sos, data_pad, axis=0)\n", + " # Cut padding to original signal\n", + " data = data_filt[npad:-npad, :]\n", + " return data, header" + ] + }, + { + "cell_type": "code", + "execution_count": 100, + "id": "4d37033d-11a3-4cfd-b0d2-62d456aed127", + "metadata": {}, + "outputs": [], + "source": [ + "def convert_2d_nifti(mask_bin, nifti_data):\n", + " # convert 4d nifti to 2d matrix using mask\n", + " nonzero_indx = np.nonzero(mask_bin)\n", + " nifti_2d = nifti_data[nonzero_indx]\n", + " return nifti_2d.T" + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "id": "30d22c16-62a9-4f5d-b0ee-4c4dd7425ae6", + "metadata": {}, + "outputs": [], + "source": [ + "def load_scans(fps, file_format, mask_fp, normalize, \n", + " bandpass, low_cut, high_cut, tr, verbose):\n", + " # scan duration in number of trs\n", + " scan_trs = []\n", + " # get # of scans\n", + " n_scans = len(fps)\n", + " # if file_format = 'nifti', load mask\n", + " if file_format == 'nifti':\n", + " mask = nb.load(mask_fp)\n", + " mask_bin = mask.get_fdata() > 0\n", + " else:\n", + " mask = None\n", + " mask_bin = None\n", + "\n", + " # initialize group matrix with zeros\n", + " group_data = initialize_matrix(fps, n_scans, file_format, \n", + " mask_bin, verbose) \n", + " print(f'loading and concatenating {n_scans} scans')\n", + " if bandpass and verbose:\n", + " print(\n", + " f'bandpass filtering of signals between {low_cut} - {high_cut} Hz '\n", + " ' will be performed'\n", + " )\n", + " # initialize counter\n", + " indx=0\n", + " # Loop through files and concatenate/append\n", + " for fp in fps:\n", + " # load file\n", + " data, header = load_file(fp, file_format, mask_bin, \n", + " bandpass, low_cut, high_cut,\n", + " tr, verbose)\n", + " # get # observations\n", + " data_n = data.shape[0]\n", + " scan_trs.append(data_n)\n", + " # Normalize data before concatenation\n", + " if normalize == 'zscore':\n", + " print(f' zscoring {fp}')\n", + " data = zscore(data, nan_policy='omit')\n", + " elif normalize == 'mean_center':\n", + " print(f' mean centering {fp}')\n", + " data = data - np.mean(data, axis=0)\n", + " # fill nans w/ 0 in regions of poor functional scan coverage\n", + " data = np.nan_to_num(data)\n", + " # fill group matrix with subject data\n", + " group_data[indx:(indx+data_n), :] = data\n", + " # increase counter by # of observations in subject scan\n", + " indx += data_n\n", + "\n", + " return group_data, mask, header, scan_trs" + ] + }, + { + "cell_type": "code", + "execution_count": 102, + "id": "590a8d95-b96c-4924-a817-26654cbec231", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "initializing matrix of size (3260, 122767)\n", + "loading and concatenating 5 scans\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010024/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010072/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010005/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010159/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-PA_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n", + " zscoring /data/SFIMJGC_Introspec/pdn/PrcsData/sub-010231/preprocessed/func/pb07_qpp/_scan_id_ses-02_task-rest_acq-AP_run-01_bold/rest2mni.b0.scale.denoise.NTRP.nii.gz\n" + ] + } + ], + "source": [ + "fps = read_input_file('/data/SFIMJGC_Introspec/cpca_vs_introspec/post_sfn_2024/config/cpca_input_output_files.txt')\n", + "#fps = read_input_file('/data/SFIMJGC_Introspec/cpca_vs_introspec/post_sfn_2024/config/cpca_input_files.txt')\n", + "\n", + "# Separate input and ouput paths\n", + "if isinstance(fps[0],list):\n", + " fps_inputs = [i[0] for i in fps]\n", + " fps_outputs = [i[1] for i in fps]\n", + "else:\n", + " fps_inputs = fps\n", + " fps_outputs = None\n", + "del fps\n", + "# master function for loading and concatenating functional scans\n", + "# parameters check\n", + "parameter_check(fps_inputs, file_format, tr, bandpass, mask_fp)\n", + "\n", + "# Pull file paths\n", + "data, mask, header, data_trs = load_scans(\n", + " fps_inputs, file_format, mask_fp, normalize, bandpass, \n", + " low_cut, high_cut, tr, verbose)" ] }, { "cell_type": "code", - "execution_count": 4, - "id": "699465fa-da41-4dac-b4a0-35db79d0422e", + "execution_count": 103, + "id": "6a282dd9-3372-4371-89b4-995487c37c04", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[652, 652, 652, 652, 652]" + ] + }, + "execution_count": 103, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "def pca(input_data, n_comps, verbose, n_iter=10):\n", - " # compute pca\n", - " print('performing PCA/CPCA')\n", - " # get number of observations\n", - " n_samples = input_data.shape[0]\n", - " print(' number of samples = %d' % n_samples)\n", - " print(' input_data.shape[1] = %d' % input_data.shape[1])\n", - " #matrix_rank = np.linalg.matrix_rank(input_data)\n", - " #print(' rank of input matrix = % s' % str(matrix_rank))\n", - " # fbpca pca\n", - " (U, s, Va) = fbpca.pca(input_data, k=n_comps, n_iter=n_iter)\n", - " # calc explained variance\n", - " explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1]\n", - " total_var = explained_variance_.sum()\n", - " # compute PC scores\n", - " pc_scores = input_data @ Va.T\n", - " # get loadings from eigenvectors\n", - " loadings = Va.T @ np.diag(s)\n", - " loadings /= np.sqrt(input_data.shape[0]-1)\n", - " # package outputs\n", - " output_dict = {'U': U,\n", - " 's': s,\n", - " 'Va': Va,\n", - " 'loadings': loadings.T,\n", - " 'exp_var': explained_variance_,\n", - " 'pc_scores': pc_scores,\n", - " 'n_samples': n_samples,\n", - " 'n_positions': input_data.shape[1],\n", - " 'total_var': total_var}\n", - " return output_dict" + "data_trs" ] }, { @@ -360,6 +1254,1009 @@ "metadata": {}, "outputs": [], "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c2b9c65-98e9-41c8-9997-77b0850ea251", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f285c94-8ce8-4af7-bcd6-440506151aa5", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7b7d14b-851d-419f-b328-1dfac9b3c78a", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e67332a-8609-479c-bc45-fa273313ece3", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd830d37-d25b-4d78-b9e6-6a68c9d104dd", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "712ce15b-b6e6-4a99-b704-65b3e8cca306", + "metadata": {}, + "outputs": [], + "source": [ + "pca_out_path = '/data/SFIMJGC_Introspec/cpca_vs_introspec/sw/complex_pca/Group_Level_133s_CPCA.on_denoise.SUBJAWARE.b0.mat'" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "752468d8-8a4e-4f4a-93c8-35929920c7ab", + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "(function(root) {\n", + " function now() {\n", + " return new Date();\n", + " }\n", + "\n", + " var force = true;\n", + " var py_version = '3.4.1'.replace('rc', '-rc.').replace('.dev', '-dev.');\n", + " var reloading = false;\n", + " var Bokeh = root.Bokeh;\n", + "\n", + " if (typeof (root._bokeh_timeout) === \"undefined\" || force) {\n", + " root._bokeh_timeout = Date.now() + 5000;\n", + " root._bokeh_failed_load = false;\n", + " }\n", + "\n", + " function run_callbacks() {\n", + " try {\n", + " root._bokeh_onload_callbacks.forEach(function(callback) {\n", + " if (callback != null)\n", + " callback();\n", + " });\n", + " } finally {\n", + " delete root._bokeh_onload_callbacks;\n", + " }\n", + " console.debug(\"Bokeh: all callbacks have finished\");\n", + " }\n", + "\n", + " function load_libs(css_urls, js_urls, js_modules, js_exports, callback) {\n", + " if (css_urls == null) css_urls = [];\n", + " if (js_urls == null) js_urls = [];\n", + " if (js_modules == null) js_modules = [];\n", + " if (js_exports == null) js_exports = {};\n", + "\n", + " root._bokeh_onload_callbacks.push(callback);\n", + "\n", + " if (root._bokeh_is_loading > 0) {\n", + " console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n", + " return null;\n", + " }\n", + " if (js_urls.length === 0 && js_modules.length === 0 && Object.keys(js_exports).length === 0) {\n", + " run_callbacks();\n", + " return null;\n", + " }\n", + " if (!reloading) {\n", + " console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n", + " }\n", + "\n", + " function on_load() {\n", + " root._bokeh_is_loading--;\n", + " if (root._bokeh_is_loading === 0) {\n", + " console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n", + " run_callbacks()\n", + " }\n", + " }\n", + " window._bokeh_on_load = on_load\n", + "\n", + " function on_error() {\n", + " console.error(\"failed to load \" + url);\n", + " }\n", + "\n", + " var skip = [];\n", + " if (window.requirejs) {\n", + " window.requirejs.config({'packages': {}, 'paths': {}, 'shim': {}});\n", + " root._bokeh_is_loading = css_urls.length + 0;\n", + " } else {\n", + " root._bokeh_is_loading = css_urls.length + js_urls.length + js_modules.length + Object.keys(js_exports).length;\n", + " }\n", + "\n", + " var existing_stylesheets = []\n", + " var links = document.getElementsByTagName('link')\n", + " for (var i = 0; i < links.length; i++) {\n", + " var link = links[i]\n", + " if (link.href != null) {\n", + "\texisting_stylesheets.push(link.href)\n", + " }\n", + " }\n", + " for (var i = 0; i < css_urls.length; i++) {\n", + " var url = css_urls[i];\n", + " if (existing_stylesheets.indexOf(url) !== -1) {\n", + "\ton_load()\n", + "\tcontinue;\n", + " }\n", + " const element = document.createElement(\"link\");\n", + " element.onload = on_load;\n", + " element.onerror = on_error;\n", + " element.rel = \"stylesheet\";\n", + " element.type = \"text/css\";\n", + " element.href = url;\n", + " console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n", + " document.body.appendChild(element);\n", + " } var existing_scripts = []\n", + " var scripts = document.getElementsByTagName('script')\n", + " for (var i = 0; i < scripts.length; i++) {\n", + " var script = scripts[i]\n", + " if (script.src != null) {\n", + "\texisting_scripts.push(script.src)\n", + " }\n", + " }\n", + " for (var i = 0; i < js_urls.length; i++) {\n", + " var url = js_urls[i];\n", + " if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n", + "\tif (!window.requirejs) {\n", + "\t on_load();\n", + "\t}\n", + "\tcontinue;\n", + " }\n", + " var element = document.createElement('script');\n", + " element.onload = on_load;\n", + " element.onerror = on_error;\n", + " element.async = false;\n", + " element.src = url;\n", + " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", + " document.head.appendChild(element);\n", + " }\n", + " for (var i = 0; i < js_modules.length; i++) {\n", + " var url = js_modules[i];\n", + " if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n", + "\tif (!window.requirejs) {\n", + "\t on_load();\n", + "\t}\n", + "\tcontinue;\n", + " }\n", + " var element = document.createElement('script');\n", + " element.onload = on_load;\n", + " element.onerror = on_error;\n", + " element.async = false;\n", + " element.src = url;\n", + " element.type = \"module\";\n", + " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", + " document.head.appendChild(element);\n", + " }\n", + " for (const name in js_exports) {\n", + " var url = js_exports[name];\n", + " if (skip.indexOf(url) >= 0 || root[name] != null) {\n", + "\tif (!window.requirejs) {\n", + "\t on_load();\n", + "\t}\n", + "\tcontinue;\n", + " }\n", + " var element = document.createElement('script');\n", + " element.onerror = on_error;\n", + " element.async = false;\n", + " element.type = \"module\";\n", + " console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", + " element.textContent = `\n", + " import ${name} from \"${url}\"\n", + " window.${name} = ${name}\n", + " window._bokeh_on_load()\n", + " `\n", + " document.head.appendChild(element);\n", + " }\n", + " if (!js_urls.length && !js_modules.length) {\n", + " on_load()\n", + " }\n", + " };\n", + "\n", + " function inject_raw_css(css) {\n", + " const element = document.createElement(\"style\");\n", + " element.appendChild(document.createTextNode(css));\n", + " document.body.appendChild(element);\n", + " }\n", + "\n", + " var js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.4.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.4.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.4.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.4.1.min.js\", \"https://cdn.holoviz.org/panel/1.4.4/dist/panel.min.js\"];\n", + " var js_modules = [];\n", + " var js_exports = {};\n", + " var css_urls = [];\n", + " var inline_js = [ function(Bokeh) {\n", + " Bokeh.set_log_level(\"info\");\n", + " },\n", + "function(Bokeh) {} // ensure no trailing comma for IE\n", + " ];\n", + "\n", + " function run_inline_js() {\n", + " if ((root.Bokeh !== undefined) || (force === true)) {\n", + " for (var i = 0; i < inline_js.length; i++) {\n", + "\ttry {\n", + " inline_js[i].call(root, root.Bokeh);\n", + "\t} catch(e) {\n", + "\t if (!reloading) {\n", + "\t throw e;\n", + "\t }\n", + "\t}\n", + " }\n", + " // Cache old bokeh versions\n", + " if (Bokeh != undefined && !reloading) {\n", + "\tvar NewBokeh = root.Bokeh;\n", + "\tif (Bokeh.versions === undefined) {\n", + "\t Bokeh.versions = new Map();\n", + "\t}\n", + "\tif (NewBokeh.version !== Bokeh.version) {\n", + "\t Bokeh.versions.set(NewBokeh.version, NewBokeh)\n", + "\t}\n", + "\troot.Bokeh = Bokeh;\n", + " }} else if (Date.now() < root._bokeh_timeout) {\n", + " setTimeout(run_inline_js, 100);\n", + " } else if (!root._bokeh_failed_load) {\n", + " console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n", + " root._bokeh_failed_load = true;\n", + " }\n", + " root._bokeh_is_initializing = false\n", + " }\n", + "\n", + " function load_or_wait() {\n", + " // Implement a backoff loop that tries to ensure we do not load multiple\n", + " // versions of Bokeh and its dependencies at the same time.\n", + " // In recent versions we use the root._bokeh_is_initializing flag\n", + " // to determine whether there is an ongoing attempt to initialize\n", + " // bokeh, however for backward compatibility we also try to ensure\n", + " // that we do not start loading a newer (Panel>=1.0 and Bokeh>3) version\n", + " // before older versions are fully initialized.\n", + " if (root._bokeh_is_initializing && Date.now() > root._bokeh_timeout) {\n", + " root._bokeh_is_initializing = false;\n", + " root._bokeh_onload_callbacks = undefined;\n", + " console.log(\"Bokeh: BokehJS was loaded multiple times but one version failed to initialize.\");\n", + " load_or_wait();\n", + " } else if (root._bokeh_is_initializing || (typeof root._bokeh_is_initializing === \"undefined\" && root._bokeh_onload_callbacks !== undefined)) {\n", + " setTimeout(load_or_wait, 100);\n", + " } else {\n", + " root._bokeh_is_initializing = true\n", + " root._bokeh_onload_callbacks = []\n", + " var bokeh_loaded = Bokeh != null && (Bokeh.version === py_version || (Bokeh.versions !== undefined && Bokeh.versions.has(py_version)));\n", + " if (!reloading && !bokeh_loaded) {\n", + "\troot.Bokeh = undefined;\n", + " }\n", + " load_libs(css_urls, js_urls, js_modules, js_exports, function() {\n", + "\tconsole.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n", + "\trun_inline_js();\n", + " });\n", + " }\n", + " }\n", + " // Give older versions of the autoload script a head-start to ensure\n", + " // they initialize before we start loading newer version.\n", + " setTimeout(load_or_wait, 100)\n", + "}(window));" + ], + "application/vnd.holoviews_load.v0+json": "(function(root) {\n function now() {\n return new Date();\n }\n\n var force = true;\n var py_version = '3.4.1'.replace('rc', '-rc.').replace('.dev', '-dev.');\n var reloading = false;\n var Bokeh = root.Bokeh;\n\n if (typeof (root._bokeh_timeout) === \"undefined\" || force) {\n root._bokeh_timeout = Date.now() + 5000;\n root._bokeh_failed_load = false;\n }\n\n function run_callbacks() {\n try {\n root._bokeh_onload_callbacks.forEach(function(callback) {\n if (callback != null)\n callback();\n });\n } finally {\n delete root._bokeh_onload_callbacks;\n }\n console.debug(\"Bokeh: all callbacks have finished\");\n }\n\n function load_libs(css_urls, js_urls, js_modules, js_exports, callback) {\n if (css_urls == null) css_urls = [];\n if (js_urls == null) js_urls = [];\n if (js_modules == null) js_modules = [];\n if (js_exports == null) js_exports = {};\n\n root._bokeh_onload_callbacks.push(callback);\n\n if (root._bokeh_is_loading > 0) {\n console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n return null;\n }\n if (js_urls.length === 0 && js_modules.length === 0 && Object.keys(js_exports).length === 0) {\n run_callbacks();\n return null;\n }\n if (!reloading) {\n console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n }\n\n function on_load() {\n root._bokeh_is_loading--;\n if (root._bokeh_is_loading === 0) {\n console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n run_callbacks()\n }\n }\n window._bokeh_on_load = on_load\n\n function on_error() {\n console.error(\"failed to load \" + url);\n }\n\n var skip = [];\n if (window.requirejs) {\n window.requirejs.config({'packages': {}, 'paths': {}, 'shim': {}});\n root._bokeh_is_loading = css_urls.length + 0;\n } else {\n root._bokeh_is_loading = css_urls.length + js_urls.length + js_modules.length + Object.keys(js_exports).length;\n }\n\n var existing_stylesheets = []\n var links = document.getElementsByTagName('link')\n for (var i = 0; i < links.length; i++) {\n var link = links[i]\n if (link.href != null) {\n\texisting_stylesheets.push(link.href)\n }\n }\n for (var i = 0; i < css_urls.length; i++) {\n var url = css_urls[i];\n if (existing_stylesheets.indexOf(url) !== -1) {\n\ton_load()\n\tcontinue;\n }\n const element = document.createElement(\"link\");\n element.onload = on_load;\n element.onerror = on_error;\n element.rel = \"stylesheet\";\n element.type = \"text/css\";\n element.href = url;\n console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n document.body.appendChild(element);\n } var existing_scripts = []\n var scripts = document.getElementsByTagName('script')\n for (var i = 0; i < scripts.length; i++) {\n var script = scripts[i]\n if (script.src != null) {\n\texisting_scripts.push(script.src)\n }\n }\n for (var i = 0; i < js_urls.length; i++) {\n var url = js_urls[i];\n if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n\tif (!window.requirejs) {\n\t on_load();\n\t}\n\tcontinue;\n }\n var element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error;\n element.async = false;\n element.src = url;\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n document.head.appendChild(element);\n }\n for (var i = 0; i < js_modules.length; i++) {\n var url = js_modules[i];\n if (skip.indexOf(url) !== -1 || existing_scripts.indexOf(url) !== -1) {\n\tif (!window.requirejs) {\n\t on_load();\n\t}\n\tcontinue;\n }\n var element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error;\n element.async = false;\n element.src = url;\n element.type = \"module\";\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n document.head.appendChild(element);\n }\n for (const name in js_exports) {\n var url = js_exports[name];\n if (skip.indexOf(url) >= 0 || root[name] != null) {\n\tif (!window.requirejs) {\n\t on_load();\n\t}\n\tcontinue;\n }\n var element = document.createElement('script');\n element.onerror = on_error;\n element.async = false;\n element.type = \"module\";\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n element.textContent = `\n import ${name} from \"${url}\"\n window.${name} = ${name}\n window._bokeh_on_load()\n `\n document.head.appendChild(element);\n }\n if (!js_urls.length && !js_modules.length) {\n on_load()\n }\n };\n\n function inject_raw_css(css) {\n const element = document.createElement(\"style\");\n element.appendChild(document.createTextNode(css));\n document.body.appendChild(element);\n }\n\n var js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.4.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.4.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.4.1.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.4.1.min.js\", \"https://cdn.holoviz.org/panel/1.4.4/dist/panel.min.js\"];\n var js_modules = [];\n var js_exports = {};\n var css_urls = [];\n var inline_js = [ function(Bokeh) {\n Bokeh.set_log_level(\"info\");\n },\nfunction(Bokeh) {} // ensure no trailing comma for IE\n ];\n\n function run_inline_js() {\n if ((root.Bokeh !== undefined) || (force === true)) {\n for (var i = 0; i < inline_js.length; i++) {\n\ttry {\n inline_js[i].call(root, root.Bokeh);\n\t} catch(e) {\n\t if (!reloading) {\n\t throw e;\n\t }\n\t}\n }\n // Cache old bokeh versions\n if (Bokeh != undefined && !reloading) {\n\tvar NewBokeh = root.Bokeh;\n\tif (Bokeh.versions === undefined) {\n\t Bokeh.versions = new Map();\n\t}\n\tif (NewBokeh.version !== Bokeh.version) {\n\t Bokeh.versions.set(NewBokeh.version, NewBokeh)\n\t}\n\troot.Bokeh = Bokeh;\n }} else if (Date.now() < root._bokeh_timeout) {\n setTimeout(run_inline_js, 100);\n } else if (!root._bokeh_failed_load) {\n console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n root._bokeh_failed_load = true;\n }\n root._bokeh_is_initializing = false\n }\n\n function load_or_wait() {\n // Implement a backoff loop that tries to ensure we do not load multiple\n // versions of Bokeh and its dependencies at the same time.\n // In recent versions we use the root._bokeh_is_initializing flag\n // to determine whether there is an ongoing attempt to initialize\n // bokeh, however for backward compatibility we also try to ensure\n // that we do not start loading a newer (Panel>=1.0 and Bokeh>3) version\n // before older versions are fully initialized.\n if (root._bokeh_is_initializing && Date.now() > root._bokeh_timeout) {\n root._bokeh_is_initializing = false;\n root._bokeh_onload_callbacks = undefined;\n console.log(\"Bokeh: BokehJS was loaded multiple times but one version failed to initialize.\");\n load_or_wait();\n } else if (root._bokeh_is_initializing || (typeof root._bokeh_is_initializing === \"undefined\" && root._bokeh_onload_callbacks !== undefined)) {\n setTimeout(load_or_wait, 100);\n } else {\n root._bokeh_is_initializing = true\n root._bokeh_onload_callbacks = []\n var bokeh_loaded = Bokeh != null && (Bokeh.version === py_version || (Bokeh.versions !== undefined && Bokeh.versions.has(py_version)));\n if (!reloading && !bokeh_loaded) {\n\troot.Bokeh = undefined;\n }\n load_libs(css_urls, js_urls, js_modules, js_exports, function() {\n\tconsole.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n\trun_inline_js();\n });\n }\n }\n // Give older versions of the autoload script a head-start to ensure\n // they initialize before we start loading newer version.\n setTimeout(load_or_wait, 100)\n}(window));" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/javascript": [ + "\n", + "if ((window.PyViz === undefined) || (window.PyViz instanceof HTMLElement)) {\n", + " window.PyViz = {comms: {}, comm_status:{}, kernels:{}, receivers: {}, plot_index: []}\n", + "}\n", + "\n", + "\n", + " function JupyterCommManager() {\n", + " }\n", + "\n", + " JupyterCommManager.prototype.register_target = function(plot_id, comm_id, msg_handler) {\n", + " if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n", + " var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n", + " comm_manager.register_target(comm_id, function(comm) {\n", + " comm.on_msg(msg_handler);\n", + " });\n", + " } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n", + " window.PyViz.kernels[plot_id].registerCommTarget(comm_id, function(comm) {\n", + " comm.onMsg = msg_handler;\n", + " });\n", + " } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n", + " google.colab.kernel.comms.registerTarget(comm_id, (comm) => {\n", + " var messages = comm.messages[Symbol.asyncIterator]();\n", + " function processIteratorResult(result) {\n", + " var message = result.value;\n", + " console.log(message)\n", + " var content = {data: message.data, comm_id};\n", + " var buffers = []\n", + " for (var buffer of message.buffers || []) {\n", + " buffers.push(new DataView(buffer))\n", + " }\n", + " var metadata = message.metadata || {};\n", + " var msg = {content, buffers, metadata}\n", + " msg_handler(msg);\n", + " return messages.next().then(processIteratorResult);\n", + " }\n", + " return messages.next().then(processIteratorResult);\n", + " })\n", + " }\n", + " }\n", + "\n", + " JupyterCommManager.prototype.get_client_comm = function(plot_id, comm_id, msg_handler) {\n", + " if (comm_id in window.PyViz.comms) {\n", + " return window.PyViz.comms[comm_id];\n", + " } else if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n", + " var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n", + " var comm = comm_manager.new_comm(comm_id, {}, {}, {}, comm_id);\n", + " if (msg_handler) {\n", + " comm.on_msg(msg_handler);\n", + " }\n", + " } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n", + " var comm = window.PyViz.kernels[plot_id].connectToComm(comm_id);\n", + " comm.open();\n", + " if (msg_handler) {\n", + " comm.onMsg = msg_handler;\n", + " }\n", + " } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n", + " var comm_promise = google.colab.kernel.comms.open(comm_id)\n", + " comm_promise.then((comm) => {\n", + " window.PyViz.comms[comm_id] = comm;\n", + " if (msg_handler) {\n", + " var messages = comm.messages[Symbol.asyncIterator]();\n", + " function processIteratorResult(result) {\n", + " var message = result.value;\n", + " var content = {data: message.data};\n", + " var metadata = message.metadata || {comm_id};\n", + " var msg = {content, metadata}\n", + " msg_handler(msg);\n", + " return messages.next().then(processIteratorResult);\n", + " }\n", + " return messages.next().then(processIteratorResult);\n", + " }\n", + " }) \n", + " var sendClosure = (data, metadata, buffers, disposeOnDone) => {\n", + " return comm_promise.then((comm) => {\n", + " comm.send(data, metadata, buffers, disposeOnDone);\n", + " });\n", + " };\n", + " var comm = {\n", + " send: sendClosure\n", + " };\n", + " }\n", + " window.PyViz.comms[comm_id] = comm;\n", + " return comm;\n", + " }\n", + " window.PyViz.comm_manager = new JupyterCommManager();\n", + " \n", + "\n", + "\n", + "var JS_MIME_TYPE = 'application/javascript';\n", + "var HTML_MIME_TYPE = 'text/html';\n", + "var EXEC_MIME_TYPE = 'application/vnd.holoviews_exec.v0+json';\n", + "var CLASS_NAME = 'output';\n", + "\n", + "/**\n", + " * Render data to the DOM node\n", + " */\n", + "function render(props, node) {\n", + " var div = document.createElement(\"div\");\n", + " var script = document.createElement(\"script\");\n", + " node.appendChild(div);\n", + " node.appendChild(script);\n", + "}\n", + "\n", + "/**\n", + " * Handle when a new output is added\n", + " */\n", + "function handle_add_output(event, handle) {\n", + " var output_area = handle.output_area;\n", + " var output = handle.output;\n", + " if ((output.data == undefined) || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n", + " return\n", + " }\n", + " var id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n", + " var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n", + " if (id !== undefined) {\n", + " var nchildren = toinsert.length;\n", + " var html_node = toinsert[nchildren-1].children[0];\n", + " html_node.innerHTML = output.data[HTML_MIME_TYPE];\n", + " var scripts = [];\n", + " var nodelist = html_node.querySelectorAll(\"script\");\n", + " for (var i in nodelist) {\n", + " if (nodelist.hasOwnProperty(i)) {\n", + " scripts.push(nodelist[i])\n", + " }\n", + " }\n", + "\n", + " scripts.forEach( function (oldScript) {\n", + " var newScript = document.createElement(\"script\");\n", + " var attrs = [];\n", + " var nodemap = oldScript.attributes;\n", + " for (var j in nodemap) {\n", + " if (nodemap.hasOwnProperty(j)) {\n", + " attrs.push(nodemap[j])\n", + " }\n", + " }\n", + " attrs.forEach(function(attr) { newScript.setAttribute(attr.name, attr.value) });\n", + " newScript.appendChild(document.createTextNode(oldScript.innerHTML));\n", + " oldScript.parentNode.replaceChild(newScript, oldScript);\n", + " });\n", + " if (JS_MIME_TYPE in output.data) {\n", + " toinsert[nchildren-1].children[1].textContent = output.data[JS_MIME_TYPE];\n", + " }\n", + " output_area._hv_plot_id = id;\n", + " if ((window.Bokeh !== undefined) && (id in Bokeh.index)) {\n", + " window.PyViz.plot_index[id] = Bokeh.index[id];\n", + " } else {\n", + " window.PyViz.plot_index[id] = null;\n", + " }\n", + " } else if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n", + " var bk_div = document.createElement(\"div\");\n", + " bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n", + " var script_attrs = bk_div.children[0].attributes;\n", + " for (var i = 0; i < script_attrs.length; i++) {\n", + " toinsert[toinsert.length - 1].childNodes[1].setAttribute(script_attrs[i].name, script_attrs[i].value);\n", + " }\n", + " // store reference to server id on output_area\n", + " output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n", + " }\n", + "}\n", + "\n", + "/**\n", + " * Handle when an output is cleared or removed\n", + " */\n", + "function handle_clear_output(event, handle) {\n", + " var id = handle.cell.output_area._hv_plot_id;\n", + " var server_id = handle.cell.output_area._bokeh_server_id;\n", + " if (((id === undefined) || !(id in PyViz.plot_index)) && (server_id !== undefined)) { return; }\n", + " var comm = window.PyViz.comm_manager.get_client_comm(\"hv-extension-comm\", \"hv-extension-comm\", function () {});\n", + " if (server_id !== null) {\n", + " comm.send({event_type: 'server_delete', 'id': server_id});\n", + " return;\n", + " } else if (comm !== null) {\n", + " comm.send({event_type: 'delete', 'id': id});\n", + " }\n", + " delete PyViz.plot_index[id];\n", + " if ((window.Bokeh !== undefined) & (id in window.Bokeh.index)) {\n", + " var doc = window.Bokeh.index[id].model.document\n", + " doc.clear();\n", + " const i = window.Bokeh.documents.indexOf(doc);\n", + " if (i > -1) {\n", + " window.Bokeh.documents.splice(i, 1);\n", + " }\n", + " }\n", + "}\n", + "\n", + "/**\n", + " * Handle kernel restart event\n", + " */\n", + "function handle_kernel_cleanup(event, handle) {\n", + " delete PyViz.comms[\"hv-extension-comm\"];\n", + " window.PyViz.plot_index = {}\n", + "}\n", + "\n", + "/**\n", + " * Handle update_display_data messages\n", + " */\n", + "function handle_update_output(event, handle) {\n", + " handle_clear_output(event, {cell: {output_area: handle.output_area}})\n", + " handle_add_output(event, handle)\n", + "}\n", + "\n", + "function register_renderer(events, OutputArea) {\n", + " function append_mime(data, metadata, element) {\n", + " // create a DOM node to render to\n", + " var toinsert = this.create_output_subarea(\n", + " metadata,\n", + " CLASS_NAME,\n", + " EXEC_MIME_TYPE\n", + " );\n", + " this.keyboard_manager.register_events(toinsert);\n", + " // Render to node\n", + " var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n", + " render(props, toinsert[0]);\n", + " element.append(toinsert);\n", + " return toinsert\n", + " }\n", + "\n", + " events.on('output_added.OutputArea', handle_add_output);\n", + " events.on('output_updated.OutputArea', handle_update_output);\n", + " events.on('clear_output.CodeCell', handle_clear_output);\n", + " events.on('delete.Cell', handle_clear_output);\n", + " events.on('kernel_ready.Kernel', handle_kernel_cleanup);\n", + "\n", + " OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n", + " safe: true,\n", + " index: 0\n", + " });\n", + "}\n", + "\n", + "if (window.Jupyter !== undefined) {\n", + " try {\n", + " var events = require('base/js/events');\n", + " var OutputArea = require('notebook/js/outputarea').OutputArea;\n", + " if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n", + " register_renderer(events, OutputArea);\n", + " }\n", + " } catch(err) {\n", + " }\n", + "}\n" + ], + "application/vnd.holoviews_load.v0+json": "\nif ((window.PyViz === undefined) || (window.PyViz instanceof HTMLElement)) {\n window.PyViz = {comms: {}, comm_status:{}, kernels:{}, receivers: {}, plot_index: []}\n}\n\n\n function JupyterCommManager() {\n }\n\n JupyterCommManager.prototype.register_target = function(plot_id, comm_id, msg_handler) {\n if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n comm_manager.register_target(comm_id, function(comm) {\n comm.on_msg(msg_handler);\n });\n } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n window.PyViz.kernels[plot_id].registerCommTarget(comm_id, function(comm) {\n comm.onMsg = msg_handler;\n });\n } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n google.colab.kernel.comms.registerTarget(comm_id, (comm) => {\n var messages = comm.messages[Symbol.asyncIterator]();\n function processIteratorResult(result) {\n var message = result.value;\n console.log(message)\n var content = {data: message.data, comm_id};\n var buffers = []\n for (var buffer of message.buffers || []) {\n buffers.push(new DataView(buffer))\n }\n var metadata = message.metadata || {};\n var msg = {content, buffers, metadata}\n msg_handler(msg);\n return messages.next().then(processIteratorResult);\n }\n return messages.next().then(processIteratorResult);\n })\n }\n }\n\n JupyterCommManager.prototype.get_client_comm = function(plot_id, comm_id, msg_handler) {\n if (comm_id in window.PyViz.comms) {\n return window.PyViz.comms[comm_id];\n } else if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n var comm = comm_manager.new_comm(comm_id, {}, {}, {}, comm_id);\n if (msg_handler) {\n comm.on_msg(msg_handler);\n }\n } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n var comm = window.PyViz.kernels[plot_id].connectToComm(comm_id);\n comm.open();\n if (msg_handler) {\n comm.onMsg = msg_handler;\n }\n } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n var comm_promise = google.colab.kernel.comms.open(comm_id)\n comm_promise.then((comm) => {\n window.PyViz.comms[comm_id] = comm;\n if (msg_handler) {\n var messages = comm.messages[Symbol.asyncIterator]();\n function processIteratorResult(result) {\n var message = result.value;\n var content = {data: message.data};\n var metadata = message.metadata || {comm_id};\n var msg = {content, metadata}\n msg_handler(msg);\n return messages.next().then(processIteratorResult);\n }\n return messages.next().then(processIteratorResult);\n }\n }) \n var sendClosure = (data, metadata, buffers, disposeOnDone) => {\n return comm_promise.then((comm) => {\n comm.send(data, metadata, buffers, disposeOnDone);\n });\n };\n var comm = {\n send: sendClosure\n };\n }\n window.PyViz.comms[comm_id] = comm;\n return comm;\n }\n window.PyViz.comm_manager = new JupyterCommManager();\n \n\n\nvar JS_MIME_TYPE = 'application/javascript';\nvar HTML_MIME_TYPE = 'text/html';\nvar EXEC_MIME_TYPE = 'application/vnd.holoviews_exec.v0+json';\nvar CLASS_NAME = 'output';\n\n/**\n * Render data to the DOM node\n */\nfunction render(props, node) {\n var div = document.createElement(\"div\");\n var script = document.createElement(\"script\");\n node.appendChild(div);\n node.appendChild(script);\n}\n\n/**\n * Handle when a new output is added\n */\nfunction handle_add_output(event, handle) {\n var output_area = handle.output_area;\n var output = handle.output;\n if ((output.data == undefined) || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n return\n }\n var id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n if (id !== undefined) {\n var nchildren = toinsert.length;\n var html_node = toinsert[nchildren-1].children[0];\n html_node.innerHTML = output.data[HTML_MIME_TYPE];\n var scripts = [];\n var nodelist = html_node.querySelectorAll(\"script\");\n for (var i in nodelist) {\n if (nodelist.hasOwnProperty(i)) {\n scripts.push(nodelist[i])\n }\n }\n\n scripts.forEach( function (oldScript) {\n var newScript = document.createElement(\"script\");\n var attrs = [];\n var nodemap = oldScript.attributes;\n for (var j in nodemap) {\n if (nodemap.hasOwnProperty(j)) {\n attrs.push(nodemap[j])\n }\n }\n attrs.forEach(function(attr) { newScript.setAttribute(attr.name, attr.value) });\n newScript.appendChild(document.createTextNode(oldScript.innerHTML));\n oldScript.parentNode.replaceChild(newScript, oldScript);\n });\n if (JS_MIME_TYPE in output.data) {\n toinsert[nchildren-1].children[1].textContent = output.data[JS_MIME_TYPE];\n }\n output_area._hv_plot_id = id;\n if ((window.Bokeh !== undefined) && (id in Bokeh.index)) {\n window.PyViz.plot_index[id] = Bokeh.index[id];\n } else {\n window.PyViz.plot_index[id] = null;\n }\n } else if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n var bk_div = document.createElement(\"div\");\n bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n var script_attrs = bk_div.children[0].attributes;\n for (var i = 0; i < script_attrs.length; i++) {\n toinsert[toinsert.length - 1].childNodes[1].setAttribute(script_attrs[i].name, script_attrs[i].value);\n }\n // store reference to server id on output_area\n output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n }\n}\n\n/**\n * Handle when an output is cleared or removed\n */\nfunction handle_clear_output(event, handle) {\n var id = handle.cell.output_area._hv_plot_id;\n var server_id = handle.cell.output_area._bokeh_server_id;\n if (((id === undefined) || !(id in PyViz.plot_index)) && (server_id !== undefined)) { return; }\n var comm = window.PyViz.comm_manager.get_client_comm(\"hv-extension-comm\", \"hv-extension-comm\", function () {});\n if (server_id !== null) {\n comm.send({event_type: 'server_delete', 'id': server_id});\n return;\n } else if (comm !== null) {\n comm.send({event_type: 'delete', 'id': id});\n }\n delete PyViz.plot_index[id];\n if ((window.Bokeh !== undefined) & (id in window.Bokeh.index)) {\n var doc = window.Bokeh.index[id].model.document\n doc.clear();\n const i = window.Bokeh.documents.indexOf(doc);\n if (i > -1) {\n window.Bokeh.documents.splice(i, 1);\n }\n }\n}\n\n/**\n * Handle kernel restart event\n */\nfunction handle_kernel_cleanup(event, handle) {\n delete PyViz.comms[\"hv-extension-comm\"];\n window.PyViz.plot_index = {}\n}\n\n/**\n * Handle update_display_data messages\n */\nfunction handle_update_output(event, handle) {\n handle_clear_output(event, {cell: {output_area: handle.output_area}})\n handle_add_output(event, handle)\n}\n\nfunction register_renderer(events, OutputArea) {\n function append_mime(data, metadata, element) {\n // create a DOM node to render to\n var toinsert = this.create_output_subarea(\n metadata,\n CLASS_NAME,\n EXEC_MIME_TYPE\n );\n this.keyboard_manager.register_events(toinsert);\n // Render to node\n var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n render(props, toinsert[0]);\n element.append(toinsert);\n return toinsert\n }\n\n events.on('output_added.OutputArea', handle_add_output);\n events.on('output_updated.OutputArea', handle_update_output);\n events.on('clear_output.CodeCell', handle_clear_output);\n events.on('delete.Cell', handle_clear_output);\n events.on('kernel_ready.Kernel', handle_kernel_cleanup);\n\n OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n safe: true,\n index: 0\n });\n}\n\nif (window.Jupyter !== undefined) {\n try {\n var events = require('base/js/events');\n var OutputArea = require('notebook/js/outputarea').OutputArea;\n if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n register_renderer(events, OutputArea);\n }\n } catch(err) {\n }\n}\n" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ] + }, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p1002" + } + }, + "output_type": "display_data" + } + ], + "source": [ + "from scipy.io import loadmat\n", + "import pandas as pd\n", + "import numpy as np\n", + "import hvplot.pandas" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "31fe8c44-0937-417a-b917-7b2b7a97fa09", + "metadata": {}, + "outputs": [], + "source": [ + "pca_out = loadmat(pca_out_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "id": "a5ae39c1-1ccb-4d82-ae7b-9fdfe4e30a3a", + "metadata": {}, + "outputs": [], + "source": [ + "out = pd.DataFrame(np.vstack([pca_out['s'], pca_out['exp_var'], pca_out['s_modified']]).T,columns=['s','exp_var','s_modified'])\n", + "out.index.name = 'cpca_id'" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "id": "11265a7c-00c0-486a-9f4f-a28c2879b2ba", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "False" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "'exp_vadr' in pca_out" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "id": "6529591e-f0bc-40c1-9d26-65c409e80910", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.4142135623730951" + ] + }, + "execution_count": 41, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.sqrt(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "id": "dd691c27-1b7c-4e4c-97f5-b74382dbdedf", + "metadata": {}, + "outputs": [], + "source": [ + "pca_out_path = '/data/SFIMJGC_Introspec/cpca_vs_introspec/spring_2024/data/cpca_javier133/Group_Level_133s_CPCA.on_denoise.SUBJAWARE.b0_exp_var.txt'\n", + "pca_out = pd.read_csv(pca_out_path)\n", + "pca_out['exp_var'] = 100 * pca_out['exp_var']" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "id": "6a499689-ddd5-4d01-bd0f-2a4ec636a3d3", + "metadata": {}, + "outputs": [ + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + ":Overlay\n", + " .Curve.Exp_var :Curve [o] (exp_var)\n", + " .Scatter.Exp_var :Scatter [index] (exp_var)" + ] + }, + "execution_count": 70, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p1443" + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "pca_out['exp_var'].hvplot('o') * pca_out['exp_var'].hvplot.scatter(s=3,c='r')" + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "id": "f1733c35-0879-4d2f-89bb-3e0fb6407ae0", + "metadata": {}, + "outputs": [ + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + ":Overlay\n", + " .Curve.Exp_var :Curve [o] (exp_var)\n", + " .Scatter.Exp_var :Scatter [index] (exp_var)" + ] + }, + "execution_count": 75, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p1775" + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "pca_out['exp_var'].cumsum().hvplot('o') * pca_out['exp_var'].cumsum().hvplot.scatter(s=3,c='r')" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "id": "00a1ce56-a949-4ac9-a780-6ab5d47187bf", + "metadata": {}, + "outputs": [], + "source": [ + "import h5py\n", + "\n", + "def load_dict_from_hdf5(filename):\n", + " data_dict = {}\n", + " with h5py.File(filename, 'r') as h5file:\n", + " for key in h5file.keys():\n", + " data_dict[key] = h5file[key][:]\n", + " return data_dict" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "db648734-1738-4ad7-8ac8-f07d3c87c3f1", + "metadata": {}, + "outputs": [], + "source": [ + "pca_out = load_dict_from_hdf5('/data/SFIMJGC_Introspec/cpca_vs_introspec/spring_2024/data/cpca_javier133/Group_Level_133s_CPCA.on_denoise.SUBJAWARE.b0.h5')" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "dda6c4e0-3fe3-40d2-a354-4f005207630c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "dict_keys(['U', 'Va', 'eigs', 'exp_var', 'loadings', 'pc_scores', 's'])" + ] + }, + "execution_count": 57, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pca_out.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5c5990f-18f8-4590-a88e-e1b593ffe8f7", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/utils/load_write.py b/utils/load_write.py index 069739e..c8f4d11 100644 --- a/utils/load_write.py +++ b/utils/load_write.py @@ -51,7 +51,7 @@ def initialize_matrix(fps, nscans, file_format, mask, verbose): elif file_format == 'txt': n_t += sum(1 for _ in open(fp)) if verbose: - print(f'initializing matrix of size ({n_t}, {n_ts})') + print(f' + [initialize_matrix]:initializing matrix of size ({n_t}, {n_ts})') # initialize matrix with zeros matrix_init = np.zeros((n_t, n_ts)) return matrix_init @@ -120,7 +120,7 @@ def load_scans(fps, file_format, mask_fp, normalize, # initialize group matrix with zeros group_data = initialize_matrix(fps, n_scans, file_format, mask_bin, verbose) - print(f'loading and concatenating {n_scans} scans') + print(f' + [load_scans]: loading and concatenating {n_scans} scans') if bandpass and verbose: print( f'bandpass filtering of signals between {low_cut} - {high_cut} Hz ' @@ -130,6 +130,7 @@ def load_scans(fps, file_format, mask_fp, normalize, indx=0 # Loop through files and concatenate/append for n,fp in enumerate(fps): + nn = n + 1 # For the output # load file data, header = load_file(fp, file_format, mask_bin, bandpass, low_cut, high_cut, @@ -139,7 +140,7 @@ def load_scans(fps, file_format, mask_fp, normalize, scan_trs.append(data_n) # Normalize data before concatenation if normalize == 'zscore': - print(f' [{n}/{n_scans}] zscoring {fp}') + print(f' [{nn}/{n_scans}] zscoring {fp}') data = zscore(data, nan_policy='omit') elif normalize == 'mean_center': print(f' mean centering {fp}') @@ -224,13 +225,13 @@ def read_input_file(input_files): # remove extra lines, if any fps = [line for line in fps if len(line)>0] # Only separate into input output if there is more than one entry per line - print(' + [read_input_file]: Number of entries in input file: %d lines' % (len(fps))) + print(' + [read_input_file]: Number of entries in input file: %d lines' % (len(fps))) if '\t' in fps[0]: fps = [item.split('\t') for item in fps] - print(' + [read_input_file]: Number of paths per input: %d paths' % (len(fps[0]))) + print(' + [read_input_file]: Number of paths per input: %d paths' % (len(fps[0]))) elif ' ' in fps[0]: fps = [item.split(' ') for item in fps] - print(' + [read_input_file]: Number of paths per input: %d paths' % (len(fps[0]))) + print(' + [read_input_file]: Number of paths per input: %d paths' % (len(fps[0]))) return fps def write_out(data, mask, header, file_format, out_prefix): From 2772bb9e7102d00a6aa8bfd5ed6404116661f935 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 17 Jan 2025 12:49:18 -0500 Subject: [PATCH 11/13] new cpca_simplified.py --- cpca_simplified.py | 447 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 447 insertions(+) create mode 100644 cpca_simplified.py diff --git a/cpca_simplified.py b/cpca_simplified.py new file mode 100644 index 0000000..f0de3f1 --- /dev/null +++ b/cpca_simplified.py @@ -0,0 +1,447 @@ +import argparse +import numpy as np +import pandas as pd +import fbpca +import warnings +import h5py +import os #CW added +from numpy.linalg import pinv +from scipy.io import savemat +from scipy.signal import hilbert +from scipy.stats import zscore, tvar +from utils.cpca_reconstruction import cpca_recon +from utils.load_write import load_data, write_out, write_modified_scans +from xmca.tools.rotation import varimax, promax + +# def save_dict_to_hdf5(filename, data_dict): +# with h5py.File(filename, 'w') as h5file: +# for key, value in data_dict.items(): +# h5file.create_dataset(key, data=value) + +def save_dict_to_hdf5(filename, data_dict): + with h5py.File(filename, 'w') as h5file: + for key, value in data_dict.items(): + if isinstance(value, np.ndarray): + # Check if the array contains objects + if value.dtype == np.dtype('O'): + # Try to convert the object array to a numeric array if possible + try: + value = np.array(value, dtype=np.float64) + except ValueError: + print(f"Skipping key {key} - cannot convert to numeric.") + continue + h5file.create_dataset(key, data=value) + else: + print(f"Skipping key {key} - not a numpy array.") + +def hilbert_transform(input_data, verbose): + # hilbert transform + input_data = hilbert(input_data, axis=0) + return input_data.conj() + +def package_parameters(n_comps, mask_fp, file_format, + pca_type, rotate, normalize, bandpass, + low_cut, high_cut, tr): + # place input parameters into dictionary to write with results + params = { + 'n_components': n_comps, + 'mask': mask_fp, + 'file_format': file_format, + 'pca_type': pca_type, + 'rotate': rotate, + 'normalize': normalize, + 'bandpass': bandpass, + 'lowcut': low_cut, + 'highcut': high_cut, + 'tr': tr + } + if mask_fp is None: + params['mask'] = '' + if rotate is None: + params['rotate'] = '' + if not bandpass: + print('bandpass not activated --> setting all bandpass related inputs to empty') + params['bandpass'] = '' + params['lowcut'] = '' + params['highcut'] = '' + if tr is None: + params['tr'] = '' + return params + + +def pca(input_data,total_var_in_func_data, n_comps, pca_type, verbose, n_iter=2): + # compute pca + print(' + [pca]: Entereing the PCA function --> performing PCA/CPCA') + # get number of observations + n_samples = input_data.shape[0] + n_vertices = input_data.shape[1] + print(' + [pca] INFO: number of samples = %d' % n_samples) + print(' + [pca] INFO: number of vertices/voxels = %d' % n_vertices) + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## print(' + [pca] INFO: Total variance in the data = %.f' % total_var_in_func_data) + # fbpca pca + (U, s, Va) = fbpca.pca(input_data, k=n_comps, n_iter=n_iter) + # calc explained variance + # 1. Get eigs + eigs = (s ** 2) / (n_samples-1) + # 2. Compute Variance Explained from eigenvalues + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## explained_variance_ = 100*np.array([eig/total_var_in_func_data for eig in eigs]) + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## total_var = explained_variance_.sum() + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## print(' + [pca] INFO: Total variance explained in PCA = %s %%' % str(total_var)) + + # OLDER CODE THAT IS INNACURATE WHEN YOU HAVE VOXELS WITH FLAT TIMESERIES + # Assumptions: + # a) We are working with complex data --> we count the number of vertices twice + # b) Input were normalized timseries --> each vertex/voxel has a variance = 1 + #if pca_type == 'complex': + # explained_variance_ = 100*np.array([eig/(n_vertices*2) for eig in eigs]) + #if pca_type == 'real': + # explained_variance_ = 100*np.array([eig/(n_vertices) for eig in eigs]) + # Original Formulation + # explained_variance_ = ((s ** 2) / (n_samples - 1)) / input_data.shape[1] + + # 3. Compute PC scores + pc_scores = input_data @ Va.T + # get loadings from eigenvectors + loadings = Va.T @ np.diag(s) + loadings /= np.sqrt(input_data.shape[0]-1) + + # 4. Package outputs + output_dict = {'U': U, + 's': s, + 'Va': Va, + 'loadings': loadings.T, + 'exp_var': np.nan, + 'eigs': eigs, + 'pc_scores': pc_scores, + 'n_samples': n_samples, + 'n_positions': input_data.shape[1], + 'total_var': np.nan} #total_var} + #'exp_var': explained_variance_, + return output_dict + + +def rotation(pca_output, data, rotation, verbose): + print(f'applying {rotation} to PCA/CPCA loadings') + # rotate PCA weights, if specified, and recompute pc scores + if rotation == 'varimax': + rotated_weights, r_mat = varimax(pca_output['loadings'].T) + pca_output['r_mat'] = r_mat + elif rotation == 'promax': + rotated_weights, r_mat, phi_mat = promax(pca_output['loadings'].T) + pca_output['r_mat'] = r_mat + pca_output['phi_mat'] = phi_mat + # https://stats.stackexchange.com/questions/59213/how-to-compute-varimax-rotated-principal-components-in-r + # recompute pc scores + projected_scores = data @ pinv(rotated_weights).T + pca_output['loadings'] = rotated_weights.T + pca_output['pc_scores'] = projected_scores + return pca_output + + +def write_results(pca_output, pca_type, mask, file_format, + header, rotate, out_prefix, n_comps_to_save): + # write out results of pca analysis + # create output name if out_prefix is None + if out_prefix is None: + if pca_type == 'complex': + out_prefix = f'cpca' + elif pca_type == f'real': + out_prefix = 'pca' + if rotate is not None: + out_prefix += f'_{rotate}' + out_prefix += '_results' + + # Write variance explained + out_df = None + if 's_modified' in pca_output: + out_df = pd.DataFrame(np.vstack([pca_output['s'], pca_output['eigs'], pca_output['exp_var'], pca_output['s_modified']]).T,columns=['s','eigs','exp_var','s_modified']) + else: + out_df = pd.DataFrame(np.vstack([pca_output['s'], pca_output['eigs'], pca_output['exp_var']]).T,columns=['s','eigs','exp_var']) + out_df.index.name = 'cpca_id' + out_df_path = f'{out_prefix}_exp_var.txt' + out_df.to_csv(out_df_path) + print(" + [write_results]: Explained Variance Dataframe written to disk [%s]" % out_df_path) + + # get loadings + loadings = pca_output['loadings'][0:n_comps_to_save:] + print('=====================') + print(loadings.shape) + print(type(pca_output)) + print('=====================') + # Write brain maps + if file_format in ('nifti', 'cifti'): + if pca_type == 'complex': + # if complex, write out real, imaginary comps, amplitude and angles + write_out( + np.abs(loadings), mask, header, + file_format, f'{out_prefix}_magnitude' + ) + write_out( + np.real(loadings), mask, header, + file_format, f'{out_prefix}_real' + ) + write_out( + np.imag(loadings), mask, header, + file_format, f'{out_prefix}_imag' + ) + write_out( + np.angle(loadings), mask, header, + file_format, f'{out_prefix}_phase' + ) + elif pca_type == 'real': + write_out( + loadings, mask, header, file_format, out_prefix + ) + # write out pca results dictionary to .mat file + print(" + Starting to save pca_output as h5 object",end=' ') + save_dict_to_hdf5(f'{out_prefix}.h5', pca_output) + #savemat(f'{out_prefix}.mat', pca_output) + print("[DONE]") + +def run_cpca(input_files, n_comps, mask_fp, file_format, out_prefix, + pca_type, rotate, recon, normalize, bandpass, + low_cut, high_cut, tr, n_bins, verbose,recon_data,n_comps_to_remove, n_comps_to_recon,save_pca_out, calc_rank): + print('++ [run_cpca]: Entering Run cpca...') + print(' + number of components to recon = %s' % str(n_comps_to_recon)) + print(' + number of components to compute = %s' % str(n_comps)) + print(' + number of components to remove = %s' % str(n_comps_to_remove)) + print(' + recon data after component removal? %s' % str(recon_data)) + print(' + recon components separately? %s' % str(recon)) + print(' + bandpass = %s' % str(bandpass)) + print(' + verbose = %s' % str(verbose)) + print(' + calculate rank of data = %s' % str(calc_rank)) + # load dataset + print(" + [run_cpca]: Loading data into memory.....") + func_data, mask, header, func_data_trs, input_paths, out_asis_paths, out_removed_paths = load_data( + input_files, file_format, mask_fp, normalize, + bandpass, low_cut, high_cut, tr, verbose + ) + + # if pca_type is complex, compute hilbert transform + if pca_type == 'complex': + print(' + [run_cpca]: Applying Hilbert Transform ...') + func_data = hilbert_transform(func_data, verbose) + + # CW: save hilbert transformed data to confirm same with new environment + #func_dict = {} + #func_dict['func_data'] = func_data + #save_dict_to_hdf5(f'{out_prefix}_func_data.h5',func_dict) + + # Estimate and write total variance in the input data (once concatenated and hilbert) + #var_in_func_data = func_data.var(axis=0) + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## var_in_func_data = tvar(func_data, axis=0, ddof=0) # CW added + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## print(' + [run_cpca]: calculated variance') + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## total_var_in_func_data = var_in_func_data.sum() + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## print(' + [run_cpca]: summed variance') + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## total_var_in_func_data_df = pd.DataFrame(var_in_func_data,columns=['Voxelwise Variance']) + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## total_var_in_func_data_df.index.name = 'voxel_ID' + + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## total_var_in_func_data_df_path = f'{out_prefix}_{pca_type}_total_var_in_func_data.txt' + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## total_var_in_func_data_df.to_csv(total_var_in_func_data_df_path) + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## print(" + [run_cpca]: Wrote variance of original data in [%s]" % total_var_in_func_data_df_path) + + # if requested, calculate the rank of the data + if calc_rank == True: + func_data_rank = np.linalg.matrix_rank(func_data) + print(' + [run_cpca]: Input data rank = %d' % func_data_rank) + + # if n_comps not provided, set it to maximum possible + if n_comps is None: + if calc_rank == True: + n_comps = func_data_rank + print(f' + [run_cpca]: Automatically setting n_comps = rank ==> n_comps = {n_comps}') + else: + n_comps = np.min(func_data.shape) + print(f' + [run_cpca]: Automatically setting n_comps = min(func_data.shape) ==> n_comps = {n_comps}') + + # compute pca + ## JAVIER REMOVED ON DEC/16/24 TO SIMPLIFY DURING TESTS ## pca_output = pca(func_data, total_var_in_func_data, n_comps, pca_type, verbose) + pca_output = pca(func_data, None, n_comps, pca_type, verbose) + + # rotate pca weights, if specified + if rotate is not None: + print(' + [run_cpca]: Applying rotation ...') + pca_output = rotation(pca_output, func_data, rotate, verbose) + + # free memory + print(' + [run_cpca]: Freeing memory by deleting the func_data variable') + del func_data + + # if cpca, and recon_data = True, create reconstructed data after removal of complex PC components + # ================================================================================================ + if recon_data & (pca_type == 'complex') & (n_comps_to_remove is not None): + if verbose: + print(f' + [run_cpca]: Performing data reconstruction after removal of {n_comps_to_remove} components') + # Nulling the first n_comps_to_remove components + pca_output['s_modified'] = pca_output['s'].copy() + pca_output['s_modified'][:n_comps_to_remove] = 0 + # Reconstructing the data (in analytical form) + func_data_modified = np.dot(pca_output['U'] * pca_output['s_modified'], pca_output['Va']) + func_data_allcomps = np.dot(pca_output['U'] * pca_output['s'], pca_output['Va']) + # Extracting the real part of the reconstructed data + func_data_modified = np.real(func_data_modified) + func_data_allcomps = np.real(func_data_allcomps) + # Write reconstructed data to disk + write_modified_scans(func_data_modified,mask,header,file_format,func_data_trs, out_removed_paths, verbose) + write_modified_scans(func_data_allcomps,mask,header,file_format,func_data_trs, out_asis_paths, verbose) + # free memory + print(' + [run_cpca]: Freeing memory by deleting the func_data_modified variable') + del func_data_modified + del func_data_allcomps + + # if cpca, and recon=True, create reconstructed time courses of complex PC + # ======================================================================== + if recon & (pca_type == 'complex'): + if verbose: + print(' + [run_cpca]: performing CPCA component time series reconstruction') + if n_comps_to_recon > 10: + warnings.warn( + 'the # of components estimated is large, CPCA reconstruction ' + 'will create a separate file for each component. This may take ' + 'a while.' + ) + del func_data # free up memory + cpca_recon(pca_output, rotate, file_format, + mask, header, out_prefix, n_bins, n_comps_to_recon) + elif recon & (pca_type == 'real'): + warnings.warn('Time series reconstruction only available for CPCA') + + # put input parameters into PCA results dicitonary + pca_output['params'] = package_parameters( + n_comps, mask_fp, file_format, pca_type, rotate, + normalize, bandpass, low_cut, high_cut, tr + ) + # write out results + if verbose: + print(' + [run_cpca]: writing out results') + if save_pca_out: + print(' + writing pca_output dictorionary to disk') + write_results(pca_output, pca_type, mask, file_format, + header, rotate, out_prefix, n_comps_to_recon) + +if __name__ == '__main__': + """Run complex or standard principal component analysis""" + parser = argparse.ArgumentParser(description='Run CPCA or PCA analysis') + parser.add_argument('-i', '--input', + help=' file path to .txt file containing the file paths ' + 'to individual fMRI scans in nifti, cifti or 2D matrices represented ' + 'in .txt. format. The 2D matrix should be observations in rows and ' + 'columns are voxel/ROI/vertices', + required=True, + type=str) + parser.add_argument('-n', '--n_comps', + help=' Number of components from PCA', + required=False, + default=None, + type=int) + parser.add_argument('-n_comps_to_recon', '--n_comps_to_recon', + help=' Number of components from PCA to reconstruct (one at a time)', + required=False, + default=None, + type=int) + parser.add_argument('-rn', '--n_comps_to_remove', + help=' Number of components to remove when reconstructing the data', + required=False, + default=None, + type=int) + parser.add_argument('-m', '--mask', + help='path to brain mask in nifti format. Only needed ' + 'if file_format="nifti"', + default=None, + required=False, + type=str) + parser.add_argument('-f', '--file_format', + help='the file format of the individual fMRI scans ' + 'specified in input', + required=False, + default='nifti', + choices=['nifti', 'cifti', 'txt'], + type=str) + parser.add_argument('-o', '--output_prefix', + help='the output file name. Default will be to save ' + 'to current working directory with standard name', + required=False, + default=None, + type=str) + parser.add_argument('-t', '--pca_type', + help='Calculate complex or real PCA', + default='complex', + choices=['real', 'complex'], + type=str) + parser.add_argument('-r', '--rotate', + help='Whether to rotate pca weights', + default=None, + required=False, + choices=['varimax', 'promax'], + type=str) + parser.add_argument('-recon', '--recon', + help='Whether to reconstruct time courses from ' + 'complex PCA', + action='store_true', + required=False) + parser.add_argument('-recon_data', '--recon_data', + help='Whether to reconstruct the data after removal of compoments', + action='store_true', + required=False) + parser.add_argument('-save_pca_out', '--save_pca_out', + help='Save the PCA dictionary', + action='store_true', + required=False) + parser.add_argument('-norm', '--normalize', + help='Type of scan normalization before group ' + 'concatenation. It is recommend to z-score', + default='zscore', + required=False, + choices=['zscore', 'mean_center'], + type=str) + parser.add_argument('-b', '--bandpass_filter', + help='Whether to bandpass filter time course w/ ' + 'a butterworth filter', + action='store_true', + required=False) + parser.add_argument('-calc_rank', '--calc_rank', + help='Whether to estimate the rank of the input data', + action='store_true', + required=False) + parser.add_argument('-f_low', '--bandpass_filter_low', + help='Low cut frequency for bandpass filter in Hz', + required=False, + default=0.01, + type=float + ) + parser.add_argument('-f_high', '--bandpass_filter_high', + help='High cut frequency for bandpass filter in Hz', + required=False, + default=0.1, + type=float + ) + parser.add_argument('-tr', '--sampling_unit', + help='The sampling unit of the signal - i.e. the TR', + required=False, + default=None, + type=float + ) + parser.add_argument('-n_bins', '--n_recon_bins', + help='Number of phase bins for reconstruction of CPCA ' + 'components. Higher number results in finer temporal ' + 'resolution', + required=False, + default=30, + type=int + ) + parser.add_argument('-v', '--verbose_off', + help='turn off printing', + action='store_false', + required=False) + + + args_dict = vars(parser.parse_args()) + os.environ["MKL_INTERFACE_LAYER"] = "ILP64" #CW added + run_cpca(args_dict['input'], args_dict['n_comps'], args_dict['mask'], + args_dict['file_format'], args_dict['output_prefix'], + args_dict['pca_type'], args_dict['rotate'], + args_dict['recon'], args_dict['normalize'], + args_dict['bandpass_filter'], args_dict['bandpass_filter_low'], + args_dict['bandpass_filter_high'], args_dict['sampling_unit'], + args_dict['n_recon_bins'], args_dict['verbose_off'],args_dict['recon_data'],args_dict['n_comps_to_remove'], + args_dict['n_comps_to_recon'],args_dict['save_pca_out'],args_dict['calc_rank']) From 4c47ea2f067d2d1a489560815843969f4018331d Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 17 Jan 2025 12:49:39 -0500 Subject: [PATCH 12/13] new cpca_simplified_remove_components --- cpca_simplified_remove_components.py | 99 ++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 cpca_simplified_remove_components.py diff --git a/cpca_simplified_remove_components.py b/cpca_simplified_remove_components.py new file mode 100644 index 0000000..1598c63 --- /dev/null +++ b/cpca_simplified_remove_components.py @@ -0,0 +1,99 @@ +import argparse +import numpy as np +import pandas as pd +import h5py +import os #CW added +from utils.load_write import load_data, write_modified_scans + +import dask.array as da +from dask.distributed import Client + +def run_recon(input_files, pca_path, n_comps_to_remove, mask_fp, normalize, verbose): + print('++ [run_recon]: Entering Run Recon...') + print(' + Number of components to remove = %s' % str(n_comps_to_remove)) + print(' + Path to mask file = %s' % mask_fp) + print(' + Path to input/output file info file = %s' % input_files) + print(' + Path to pre-computed CPCA object = %s' % pca_path) + print(' + Verbose? = %s' % str(verbose)) + print(' + Normalize = %s' % str(normalize)) + # load dataset + print(" + [run_recon]: Loading data into memory.....") + func_data, mask, header, func_data_trs, _, _, out_removed_paths = load_data(input_files, 'nifti', mask_fp, normalize, False, None, None, None, verbose) + + # Load CPCA objects + print(" + [run_recon]: Load CPCA objects into memory with Dask....") + dask_client = Client() + print(dask_client) + + f = h5py.File(pca_path, "r") + pca_output = {key: da.from_array(f[key], chunks="auto") for key in f.keys()} + + print(' + Loading s...', end='') + s = pca_output['s'].compute() + print(' [DONE] s.shape=%s' % str(s.shape)) + + print(' + Loading Va...', end='') + Va = pca_output['Va'].compute() + print(' [DONE] Va.shape=%s' % str(Va.shape)) + + print(' + Loading U...', end='') + U = pca_output['U'].compute() + print(' [DONE] U.shape=%s' % str(U.shape)) + + # Set to Zero the components to be removed + print('++ [run_recon]: Setting to zero the components to be removed') + print(' + s = %s' % str(s)) + s_modified = s.copy() + s_modified[:n_comps_to_remove] = 0 + print(' + s_modified = %s' % str(s_modified)) + + # Reconstructing the data (in analytical form) + print('++ [run_recon]: Computing the modified data...') + func_data_modified = np.dot(U * s_modified, Va) + # Extracting the real part of the reconstructed data + print('++ [run_recon]: Extracting the real part of the modified data...') + func_data_modified = np.real(func_data_modified) + # Write reconstructed data to disk + print('++ [run_recon]: Writing modified scans to disk...') + write_modified_scans(func_data_modified,mask,header,'nifti',func_data_trs, out_removed_paths, verbose) + print('++ [run_recon]: Program ends') + +if __name__ == '__main__': + """Reconstruct data using pre-existing CPCA decomposition after removal of a given number of components""" + parser = argparse.ArgumentParser(description='Data recon from pre-computed CPCA') + parser.add_argument('-i', '--input', + help=' file path to .txt file containing the file paths ' + 'to individual fMRI scans in nifti, cifti or 2D matrices represented ' + 'in .txt. format. The 2D matrix should be observations in rows and ' + 'columns are voxel/ROI/vertices', + required=True, + type=str) + parser.add_argument('-pca','--pca_path', + help=' path to previously computed CPCA', + required=True, + type=str) + parser.add_argument('-rn', '--n_comps_to_remove', + help=' Number of components to remove when reconstructing the data', + required=False, + default=None, + type=int) + parser.add_argument('-m', '--mask', + help='path to brain mask in nifti format. Only needed ' + 'if file_format="nifti"', + default=None, + required=False, + type=str) + parser.add_argument('-v', '--verbose_off', + help='turn off printing', + action='store_false', + required=False) + parser.add_argument('-norm', '--normalize', + help='Type of scan normalization before group ' + 'concatenation. It is recommend to z-score', + default='zscore', + required=False, + choices=['zscore', 'mean_center'], + type=str) + args_dict = vars(parser.parse_args()) + os.environ["MKL_INTERFACE_LAYER"] = "ILP64" #CW added + run_recon(args_dict['input'], args_dict['pca_path'], args_dict['n_comps_to_remove'], args_dict['mask'],args_dict['normalize'],args_dict['verbose_off']) From c2f4e33ac3d112f47da49c846851d5d39b7699aa Mon Sep 17 00:00:00 2001 From: Javier Gonzalez-Castillo Date: Fri, 17 Jan 2025 12:50:04 -0500 Subject: [PATCH 13/13] new cpca_zscore.py --- cpca_zscore.py | 53 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 cpca_zscore.py diff --git a/cpca_zscore.py b/cpca_zscore.py new file mode 100644 index 0000000..7729ffe --- /dev/null +++ b/cpca_zscore.py @@ -0,0 +1,53 @@ +import argparse +import nibabel as nb +import numpy as np +from scipy.stats import zscore +from utils.load_write import load_file, write_out +import os #CW added + + +def run_zscore(input_path, output_path, mask_path): + print('++ [run_zscore]: Entering Run Recon...') + print(' + Path to mask file = %s' % mask_path) + print(' + Path to input file = %s' % input_path) + print(' + Path to output file = %s' % output_path) + # Load Mask + print('++ [run_zscore]: Load mask...') + mask = nb.load(mask_path) + mask_bin = mask.get_fdata() > 0 + + # Load input data + print('++ [run_zscore]: Load input data...') + data,header = load_file(input_path,'nifti',mask_bin,False,None,None,None,True) + data_n = data.shape[0] + + # Z-score the data across the time dimension + print('++ [run_zscore]: Z-score data across the time dimension...') + data = zscore(data, nan_policy='omit') + data = np.nan_to_num(data) + + # Write output dataset + print('++ [run_zscore]: Write output dataset...') + write_out(data,mask,header,'nifti',output_path) + + print('++ [run_zscore]: Program ends') + +if __name__ == '__main__': + """Create Z-scored version of input data (micmiking CPCA program)""" + parser = argparse.ArgumentParser(description='Create Z-scored version of the data') + parser.add_argument('-i', '--input', + help=' path to input nifti file', + required=True, + type=str) + parser.add_argument('-o','--output', + help=' path to output nifti file', + required=True, + type=str) + parser.add_argument('-m', '--mask', + help='path to mask in nifti format', + default=None, + required=True, + type=str) + args_dict = vars(parser.parse_args()) + os.environ["MKL_INTERFACE_LAYER"] = "ILP64" #CW added + run_zscore(args_dict['input'], args_dict['output'], args_dict['mask'])