Add RapidPT to SnPM

RapidPT_min/AddPaths.m

@@ -0,0 +1,11 @@
+function [] = AddPaths(RapidPTPath)
+%AddPaths Add all paths needed
+
+    % Add path mex path and grasta paths
+    PATH = RapidPTPath;
+    addpath(PATH);
+    addpath(strcat(PATH, '/util'));
+    addpath(strcat(PATH, '/include/grasta.1.2.0'));
+
+end
+

RapidPT_min/CheckSamplingRate.m

@@ -0,0 +1,19 @@
+function [ subV ] = CheckSamplingRate( N, V, sub)
+% min_subV Validate that we are sampling enough entries from the matrix
+%   * Condition: 
+%`      We want V*subV > r*log(V) 
+
+minSampling = round(N * log(V));
+subV = round(sub * V); 
+if(minSampling > subV)
+    warning('Warning: The input sub-sampling rate is too low! RapidPT requires that V*subV > N*log(V) to be able to accurately recover the Maxnull distribution. See paper reference (Sec 3.1) for more information on the theory behind this requirement.');
+    warning('We are automatically changing the sub-sampling rate such that V*sub > 2*N*log(V).');
+    minSub = 2*N*log(V)/V;
+    subV = round(minSub*V);
+    fprintf( 'Old sub-sampling rate, sub = %d \n', sub);
+    fprintf( 'New sub-sampling rate, sub = %d, subV = %d \n', minSub, subV);
+end
+
+
+end
+

RapidPT_min/LICENSE.txt

@@ -0,0 +1,13 @@
+RapidPT is provided without any warranty of fitness for any purpose. 
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.

RapidPT_min/RapidPT.m

@@ -0,0 +1,210 @@
+%% Efficient permutation testing using Matrix completion
+% % the following function computes the max and min null statistic distribution 
+% % in its current format, the code only uses t-statistics
+
+%%% Original Library:
+% % RapidPT: https://github.com/felipegb94/RapidPT
+
+%%% Corresponding papers :
+% % Accelerating Permutation Testing in Neuroimaging through Subspace Tracking: A new plugin for SnPM.
+% % F. Gutierrez-Barragan, V.K. Ithapu, C. Hinrichs, C. Maumet, S.C. Johnson, T.E. Nichols, V. Singh.
+% % Neuorimage, 2017.
+% % Speeding up Permutation Testing in Neuroimaging 
+% % C Hinrichs, VK Ithapu, Q Sun, SC Johnson, V Singh
+% % NIPS 2013
+
+%%% Arguments    
+% % 
+% %     %%% INPUTS
+% %     A structure filed with following arguments
+% %     inputs.data    :       path to mat file containing the data matrix 
+% %                                (REQUIRED) Two fields : data and labeling
+% %                                data - a matrix of size N X V
+% %                                labels - a vector of length N (2 groups)
+% %                                CONTENTS SHOULD BE NAMED "data" and "labels"                                
+% %                                (N : number of instances , V : data dimension) 
+% %     inputs.sub         :       sub-sampling rate (0 to 1) (DEFAULT = 0.05)
+% %     inputs.T           :       number of permutations (DEFAULT = 10^4)
+% %     inputs.maxrank     :       rank for estimating the low rank subspace (DEFAULT = N)
+% %     inputs.traintime   :       number of permutations for training (DEFAULT = 100)
+% %     inputs.maxCycles   :       number of cycles for training (DEFAULT = 3)
+% %     inputs.iter        :       number of iterations for matrix completion (DEFAULT = 30)  
+% %     inputs.writing     :       if 0 - outputs only maxnull (SEE BELOW) 
+% %                                if 1 - outputs maxnull, U and W (DEFAULT = 0)
+% %     inputs.save        :       path to save the outputs (DEFAULT : working folder)  
+% % 
+% %     %%% OUTPUTS
+% %     outputs.maxnull     :       estimated distribution of max Null statistic
+% %     outputs.U           :       orthogonal matrix spanning low rank subspace
+% %                                 (dimension : V X maxrank) 
+% %                                 optional output (DEFAULT : No)
+
+%%% Support codes for matrix completion
+% % GRASTA : https://sites.google.com/site/hejunzz/grasta
+% % Codes already included in the package
+
+%%% Usage 
+% %     See https://github.com/felipegb94/RapidPT
+
+function [ outputs, timings ] = RapidPT( inputs, rapidPTLibraryPath, T0 )
+
+    fprintf('Starting RapidPT...\n');
+    tTotal = tic;
+    fprintf('Adding Paths...\n');
+    AddPaths(rapidPTLibraryPath);
+    fprintf('\nStarting Preprocessing...\n');
+    fprintf('Validate Required Inputs...\n');
+    ValidateInputs(inputs);
+    data = inputs.data;
+    dataSquared = data.*data;
+    %labels = inputs.labels;
+    nGroup1 = inputs.nGroup1;
+
+
+
+    fprintf('Processing Input Parameters...\n');
+    N = size(data,1); % N: number of instances/subjects (rows in data matrix)
+    V = size(data,2); % V: Number of statistics/voxel measurements (cols) 
+%     uniqueLabels = unique(labels); % Unihttps://github.com/felipegb94/RapidPermTest.gitque labels
+%     nGroup1 = length(find(labels==uniqueLabels(1))); % Number of patients in group 1
+    nGroup2 = N - nGroup1;
+    [sub, numPermutations, maxRank, trainNum, maxCycles, iter, write] = ProcessInput(inputs, N);
+
+    fprintf('Initializing matrix completion parameters (GRASTA parameters) \n');
+    [ options, opts, opts2, status ] = InitGrastaParams(maxRank, iter, V);
+
+    fprintf('Initializing permutation matrices... \n');
+    % indexMatrix is what indexMatrix used to be..
+    [~, permutationMatrix1, permutationMatrix2] = TwoSampleGetPermutationMatrices(numPermutations, N, nGroup1);
+
+    binRes = 0.05; 
+    maxnullBins = -9:binRes:9; %% bin resolution in maxnull histogram computation
+    subV = CheckSamplingRate(N, V, sub); %% number of samples used per permutation
+%     subV = round(V*sub);
+    maxTStatistics = zeros(1, numPermutations); %% estimated max statistics for all permutations
+    minTStatistics = zeros(1, numPermutations); %% estimated max statistics for all permutations
+
+%% Training for low rank subsapace and residual priors
+
+    fprintf('\nStarting RapidPT core...\n');
+    fprintf('Training for low rank subspace and residual priors \n');
+
+    clear inputs; % Clear some memory
+    tTraining = tic;
+
+    permutationMatrix1Current = permutationMatrix1(1:trainNum,:);
+    permutationMatrix2Current = permutationMatrix2(1:trainNum,:);
+    % Calculate some full permutations for training U (A good number would be the number of labels).
+    [TCurrent] = TwoSamplePermTest(data, dataSquared, permutationMatrix1Current, permutationMatrix2Current, nGroup1, nGroup2);
+
+    framesOrder = zeros(trainNum, maxCycles);
+    for m = 1:1:maxCycles
+        framesOrder(:,m) = randperm(trainNum);
+    end
+
+    % Estimate U using subsample matrix completion methods
+    UHat = orth(randn(V,options.RANK)); 
+
+    for m = 1:1:maxCycles
+        for f = 1:1:trainNum
+            inds = randperm(V,subV); 
+            I_inds = TCurrent(framesOrder(f,m),inds)';
+            [UHat, status, opts] = grasta_stream(I_inds, inds, UHat, status, options, opts);
+            fprintf('Subspace estimation on %s cycle with %s frame \n',num2str(m),num2str(f));
+        end
+    end 
+
+    diffForNormal = zeros(trainNum,maxCycles);
+
+    for m = 1:1:maxCycles
+        Ts_ac = zeros(V,trainNum); 
+        Ts_tr = zeros(V,trainNum);
+        for f = 1:1:trainNum
+            Ts_ac(:,f) = TCurrent(framesOrder(f,m),:)';
+            inds = randperm(V,subV); 
+
+            I_inds = Ts_ac(inds,f);
+            % Time srp function
+            [s, w, ~] = admm_srp(UHat(inds,:), I_inds, opts2); 
+            sall = zeros(V,1); 
+            sall(inds) = s; 
+            Ts_tr(:,f) = (UHat*w + sall)';
+            fprintf('Training done on %s cycle with %s frame \n',num2str(m),num2str(f));
+        end
+        max_Ts_ac = max(Ts_ac,[],1); 
+        max_Ts_tr = max(Ts_tr,[],1);
+        diffForNormal(:,m) = max_Ts_ac - max_Ts_tr;
+    end
+
+    [muFit,~] = normfit(diffForNormal(:));
+
+    tTraining = toc(tTraining);
+    timings.tTraining = tTraining;
+
+%% Max Memory Check + Max Number of parallel workers
+    clear TCurrent; clear Ts_ac; clear Ts_tr;
+    c = parcluster('local'); 
+    numCoresAvail = c.NumWorkers;
+    maxBytes = 8*1024*1024*1024; % let maxMemory be 8GB
+    % Get variables with large memory footprint
+    UHatInfo = whos('UHat'); dataInfo = whos('data'); permMatrixInfo = whos('permutationMatrix1'); 
+    bytesPerWorker = UHatInfo.bytes + 2*dataInfo.bytes + 2*permMatrixInfo.bytes; % There is 2 data matrices and 2 perm matrices
+
+    maxNumWorkers = floor(maxBytes / bytesPerWorker); % If zero parfor will run serially
+    if maxNumWorkers == 1
+        maxNumWorkers = 0;
+    end
+    numWorkers = min(maxNumWorkers,numCoresAvail);
+
+%% Recovery : Filling in W and residuals for all numPermutations
+    fprintf('\n Recovering the subspace coefficients and residuals for all permutations \n');
+    tRecovery = tic;  
+    nPtmp = ones(V,1); T0 = T0'; % Initialize to ones because the first statistic is equal to T0
+    % If numWorkers==1 the following loop will simply run serially
+    parfor (i = 1:numPermutations, numWorkers)
+        inds = randperm(V,subV)'; 
+        [TCurrent] = TwoSamplePermTest(data(:,inds),...
+                                       dataSquared(:,inds),...
+                                       permutationMatrix1(i,:),...
+                                       permutationMatrix2(i,:),...
+                                       nGroup1,...
+                                       nGroup2);
+        U_inds = UHat(inds,:);  
+        [s, w, ~] = admm_srp(U_inds, TCurrent', opts2);
+        %W{i,1} = w; 
+        TRec = UHat*w;
+        TRec(inds) = TRec(inds) + s;
+        maxTStatistics(1,i) = max(TRec) + muFit;
+        minTStatistics(1,i) = min(TRec) - muFit;
+        nPtmp = nPtmp + ((TRec)>=T0); % Used by SnPM to write lP images
+        fprintf('Completion done on trial %d/%d  \n',i,numPermutations);  
+    end
+
+    % Save timings
+    timings.tRecovery = toc(tRecovery);
+
+    outputs.MaxNull = gen_hist(maxTStatistics,maxnullBins); 
+    outputs.MaxT = maxTStatistics;
+    outputs.MinT = minTStatistics;
+    outputs.nP = nPtmp';
+
+%     if write == 1 
+%         outputs.U = UHat; 
+%         %outputs.W = W; 
+%     end
+
+    fprintf('TwoSampleRapidPT Done...\n');
+    timings.tTotal = toc(tTotal);
+
+end
+
+function [tStatMatrix] = TwoSamplePermTest(data, dataSquared, permutationMatrix1, permutationMatrix2, nGroup1, nGroup2)
+
+    g1Mean = (permutationMatrix1 * data)/nGroup1;
+    g2Mean = (permutationMatrix2 * data)/nGroup2;
+    g1Var = (permutationMatrix1 * dataSquared)/(nGroup1) - (g1Mean.*g1Mean);
+    g2Var = (permutationMatrix2 * dataSquared)/(nGroup2) - (g2Mean.*g2Mean);
+    tStatMatrix = (g1Mean - g2Mean) ./ (sqrt((g1Var./(nGroup1-1)) + (g2Var./(nGroup2-1))));
+
+end
+

RapidPT_min/TwoSampleGetPermutationMatrices.m

@@ -0,0 +1,36 @@
+function [ indexMatrix, permutationMatrix1, permutationMatrix2 ] = TwoSampleGetPermutationMatrices(numPermutations, N, nGroup1 )
+%GetPermutationMatrices 
+%   * indexMatrix: Each column contains the indeces that will be used to
+%   make a permutation. Rows 1-nGroup1 are the indeces for group 1 and
+%   rows nGroup1-size(labels,1) are the indeces for group 2. This matrix is
+%   used when doing serial permutation testing where each permutation is
+%   calculated one by one.
+%
+%   * permutationMatrix1: Matrix composed of 1's and 0's. At each row the
+%   columns that contain 1's refer to the corresponding subject in the data
+%   matrix at that index and this subject will be in group 1 for that
+%   permutations.
+%
+%   * permutationMatrix2: Same as permutationMatrix2 but for group 2
+%   instead of 1.
+    rng('default');
+    rng('shuffle');
+    indexMatrix = zeros(numPermutations, N); 
+    permutationMatrix1 = zeros(numPermutations, N);
+    permutationMatrix2 = zeros(numPermutations, N);
+
+    for t = 1:1:numPermutations
+        currLabels = randperm(N);
+        currLabels1 = currLabels(1:(nGroup1));
+        currLabels2 = currLabels(nGroup1+1:end);
+        permutationMatrix1(t,currLabels1) = 1;
+        permutationMatrix2(t,currLabels2) = 1;
+        indexMatrix(t,:) = currLabels;
+    end
+
+end
+
+
+
+
+

RapidPT_min/TwoSampleRapidPT.m

@@ -0,0 +1,97 @@
+%% Efficient permutation testing using Matrix completion
+% % the following function computes the max Null statistic distribution 
+% % in its current format, the code only uses t-statistics
+
+%%% Original Library:
+% % RapidPT: https://github.com/felipegb94/RapidPT
+
+%%% Corresponding papers :
+% % Accelerating Permutation Testing in Neuroimaging through Subspace Tracking: A new plugin for SnPM.
+% % F. Gutierrez-Barragan, V.K. Ithapu, C. Hinrichs, C. Maumet, S.C. Johnson, T.E. Nichols, V. Singh.
+% % Neuorimage, 2017.
+% % Speeding up Permutation Testing in Neuroimaging 
+% % C Hinrichs, VK Ithapu, Q Sun, SC Johnson, V Singh
+% % NIPS 2013
+
+%%% Inputs: 
+%       - Data: NxV matrix - N: number of subject, V: number of voxels per subject.
+%       - numPermutations: Number of permutation to perform
+%       - nGroup1: Number of subjects in group 1
+%       - writingVal: Can be 1 or 0. 0 means that we only save the maxnull
+%         distribution produced by the algorithm. 1 means that we save the maxnull
+%         distribution AND the low-rank basis matrix and coefficient matrix with
+%         which you can recover the Permutation testing matrix
+%%% Outputs
+%       - outputs.MaxT - Maximum t-statistics at each permutation.
+%       - outputs.MaxNull - Binned MaxT statistics.
+%       - outputs.tThreshold - Calculated tThreshold for the input p-value.
+%       - timings.tTraining - Time it took for the first part of RapidPT
+%       procedure
+%       - timings.tRecovery - Time it took to recover the whole Permutation
+%       matrix.
+%       - timings.tTotal - Total time it took.
+
+function [ outputs, timings ] = TwoSampleRapidPT(Data, numPermutations, nGroup1, writingVal, RapidPTLibraryPath,T0)
+% RapidPT 
+%   Modified permutation testing algorithm described in the following paper
+%   Speeding up Permutation Testing in Neuroimaging  C Hinrichs, VK Ithapu, Q Sun, SC Johnson, V Singh, NIPS 2013
+
+    % N subjects, V voxels (or statistics)
+    [N,V] = size(Data);
+
+    % Set keys for input struct
+    dataKey = 'data';
+    %labelsKey = 'labels';
+    nGroup1Key = 'nGroup1';
+    %nGroup2Key = 'nGroup2';
+    subKey = 'sub';
+    TKey = 'T';
+    maxRankKey = 'maxRank';
+    trainNumKey = 'trainNum';
+    maxCyclesKey = 'maxCycles';
+    iterKey = 'iter';
+    writingKey = 'writing';
+
+%% Fixed values for RapidPT
+% The following parameters have been extensively tested and have produced
+% good results, independent of the dataset size and number of permutations
+
+    % Set the sub-sampling rate to low values. Inside the RapidPT function
+    % these values will be checked and IF they are too low they will be set
+    % to the correct value.
+    if(numPermutations < 10000)      
+        subVal = {0.005};
+    else
+        subVal = {0.0035};
+    end
+
+    maxCyclesVal = {3}; % Number of cycles for training.
+    iterVal = {30}; % Number of iterations for matrix completion.
+
+%% Variable inputs for RapidPT
+% These inputs depend on the dataset size.
+
+    nGroup1Val = nGroup1; % Size of group 1
+    assert(N > nGroup1, 'nGroup1 cannot be larger than the number of subjects in the data');
+    TVal = {numPermutations}; % Number of Permutations.
+    maxRankVal = {N}; % Rank for estimating the low rank subspace
+    trainNumVal = N; % Default number of training samples
+
+
+    inputs = struct(dataKey, Data,...
+                    nGroup1Key, nGroup1Val,...
+                    subKey, subVal,...
+                    TKey, TVal,...
+                    maxRankKey, maxRankVal,...
+                    trainNumKey, trainNumVal,...
+                    maxCyclesKey, maxCyclesVal,...
+                    iterKey, iterVal,...
+                    writingKey, writingVal);
+
+    [outputs, timings] = RapidPT(inputs, RapidPTLibraryPath, T0);
+
+
+end
+
+
+