diff --git a/Matlab_Scripts/Functions/RunPLS.m b/Matlab_Scripts/Functions/RunPLS.m
index 088f098..c6ab7bc 100644
--- a/Matlab_Scripts/Functions/RunPLS.m
+++ b/Matlab_Scripts/Functions/RunPLS.m
@@ -1,141 +1,141 @@
 % Both entries should have size n_subjects x n_dims WITHOUT NAN VALUES
-function [U,V,EV,sprob,BS_U,BS_V,Ub_vect,Vb_vect] = RunPLS(Functional,Behavioral)
+function [U,V,EV,sprob,BS_U,BS_V,Ub_vect,Vb_vect,permsamp,S] = RunPLS(Functional,Behavioral)
     
     Y0 = Functional;
     X0 = Behavioral;
 
     % Number of subjects in the PLS
     nSub = size(Y0,1);
     
     % normalization of X and Y
     Y=myPLS_norm(Y0,1,ones(1,nSub),2);
     X=myPLS_norm(X0,1,ones(1,nSub),2);
 
     %%% COMPUTATIONS OF ACTUAL PLS
 
     % Covariance matrix computation
     R=myPLS_cov(X,Y,1,ones(1,nSub));
     
     % SVD on the covariance matrix
     [U,S,V]=svd(R,'econ');
 
     EV = (diag(S).^2) / sum(diag(S.^2));
     
     % Number of latent variables
     NUM_LVs=min(size(S));
     
     % ICA convention: turn LVs such that max is positive
     for iLV=1:NUM_LVs
         [~,maxID] = max(abs(U(:,iLV)));
         if U(maxID,iLV)<0
             U(:,iLV)=-U(:,iLV);
             V(:,iLV)=-V(:,iLV);
         end
     end
     
     
     %%% PERMUTATION-BASED ASSESSMENT OF SIGNIFICANCE
 
     for iter=1:1000
 
         fprintf('%d ',iter);
 
         if ~mod(iter,20) 
             fprintf('\n'); 
         end
 
         % How to randomize the data
         perm_order=randperm(size(X,1));
 
         % Permuted behavioral data (untouched)
         Xp=X;
 
         % Permuted functional data (shuffled)
         Yp=Y0(perm_order,:);
 
         % Renormalization
         Yp=myPLS_norm(Yp,1,ones(1,nSub),2);
 
         % "Null" correlation matrix and SVD
         Rp=myPLS_cov(Xp,Yp,1,ones(1,nSub));
         [Up,Sp,Vp]=svd(Rp,'econ');
 
         % Procrustes transform to fit the new U (permuted) to the old one, as
         % well as the singular values
-        rotatemat = rri_bootprocrust(U, Up);
-        Up = Up * Sp * rotatemat;
-        Sp = sqrt(sum(Up.^2));
+        %rotatemat = rri_bootprocrust(U, Up);
+        %Up = Up * Sp * rotatemat;
+        %Sp = sqrt(sum(Up.^2));
 
         % Current null singular values are stored
-        permsamp(:,iter)=Sp';
+        permsamp(:,iter)=diag(Sp);
 
         % sp contains whether the current singular value exceeds the actual one
         % (real one)
         if iter==1
-            sp= (Sp'>=diag(S));
+            sp= (diag(Sp)>=diag(S));
         else
-            sp = sp + (Sp'>=diag(S));
+            sp = sp + (diag(Sp)>=diag(S));
         end
     end
 
     fprintf('\n');
 
     % The larger the sprob entries, the larger the p-value
     sprob = sp ./ (1000 + 1);
 
     for iter=1:NUM_LVs
         str{iter}=sprintf('%d (p=%5.4f)',iter,sprob(iter));
         fprintf('LV%s\n',str{iter});
     end
 
 
     %%% BOOTSTRAPPING ASSESSMENT OF STABILITY
     
     % Computes the indices to use for bootstrapping (with replacement)
     [boot_order]=rri_boot_order(nSub,1,1000);
 
     % Runs across bootstrapping folds
     for iter=1:1000
 
         fprintf('%d ',iter);
         if ~mod(iter,20)
             fprintf('\n');
         end
 
         % Takes bootstrapped behavioral data and normalizes it
         Xb=X0(boot_order(:,iter),:);
         Xb=myPLS_norm(Xb,1,ones(1,nSub),2);
 
         % Same for functional data
         Yb=Y0(boot_order(:,iter),:);
         Yb=myPLS_norm(Yb,1,ones(1,nSub),2);
 
         % Accordingly, computes R and the SVD
         Rb=myPLS_cov(Xb,Yb,1,ones(1,nSub));
         [Ub,Sb,Vb]=svd(Rb,'econ');
 
         % Procrustes transform as before to generate fitting U and V
         % bootstrapped variants
         rotatemat = rri_bootprocrust(U, Ub);
         rotatemat2 = rri_bootprocrust(V, Vb);
 
         % Full rotation
         Rot_full = (rotatemat + rotatemat2)/2;
 
         Vb = Vb * Rot_full;
         Ub = Ub * Rot_full;
 
         % Saves them all
         Ub_vect(iter,:,:) = Ub;
         Vb_vect(iter,:,:) = Vb;
     end
     
     % Generation of boostrap scores
     mu = squeeze(mean(Ub_vect,1));
     stan = squeeze(std(Ub_vect,1));
     BS_U = mu./stan;
     
     mu2 = squeeze(mean(Vb_vect,1));
     stan2 = squeeze(std(Vb_vect,1));
     BS_V = mu2./stan2;
 end
\ No newline at end of file