a59: improved lowrank and other fixes for -npccon.

palm.m

@@ -77,4 +77,4 @@ function palm(varargin)
 end
 
 % Now run what matters
-palm_backend(varargin{:});
+palm_core(varargin{:});

palm_clusterd.m

@@ -42,7 +42,7 @@
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
 % Choose an appropriate mask struct.
-if opts.NPC || opts.MV,
+if opts.npcmod || opts.MV,
     S = plm.maskinter;
 else
     if plm.nmasks == 1,

palm_clustere.m

@@ -42,7 +42,7 @@
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
 % Choose an appropriate mask struct.
-if opts.NPC || opts.MV,
+if opts.npcmod || opts.MV,
     S = plm.maskinter;
 else
     if plm.nmasks == 1,

palm_clusterm.m

@@ -42,7 +42,7 @@
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
 % Choose an appropriate mask struct.
-if opts.NPC || opts.MV,
+if opts.npcmod || opts.MV,
     S = plm.maskinter;
 else
     if plm.nmasks == 1,

palm_clustert.m

@@ -42,7 +42,7 @@
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
 % Choose an appropriate mask struct.
-if opts.NPC || opts.MV,
+if opts.npcmod || opts.MV,
     S = plm.maskinter;
 else
     if plm.nmasks == 1,

palm_csort.m → palm_competitive.m

@@ -1,12 +1,12 @@
-function [unsrtR,S,srtR] = palm_csort(X,ord,mod)
+function [unsrtR,S,srtR] = palm_competitive(X,ord,mod)
 % Sort a set of values and return their competition
 % ranks, i.e., 1224, or the modified competition ranks,
 % i.e. 1334. This makes difference only when there are
 % ties in the data. The function returns the ranks in
 % their original order as well as sorted.
 % 
 % Usage:
-% [unsrtR,S,srtR] = palm_csort(X,ord,mod)
+% [unsrtR,S,srtR] = palm_competitive(X,ord,mod)
 % 
 % Inputs:
 % - X      : 2D array with the original data. The
@@ -26,9 +26,9 @@
 %
 % Examples:
 % - To obtain the empirical cdf of a dataset in X, use:
-%   cdf   = palm_csort(X,'ascend',true)/size(X,1);
+%   cdf   = palm_competitive(X,'ascend',true)/size(X,1);
 % - To obtain the empirical p-values for each value in X, use:
-%   pvals = palm_csort(X,'descend',true)/size(X,1);
+%   pvals = palm_competitive(X,'descend',true)/size(X,1);
 % 
 % _____________________________________
 % Anderson M. Winkler
@@ -81,7 +81,7 @@
 
 % Unmodified competition ranking
 [nR,nC] = size(X);
-unsrtR = single(zeros(size(X)));
+unsrtR  = single(zeros(size(X)));
 [S,tmp] = sort(X,ord);
 [~,rev] = sort(tmp);
 srtR = repmat((1:nR)',[1 nC]);

palm_datapval.m

@@ -50,7 +50,7 @@
 if rev, % if small G are significant
 
     % Sort the data and compute the empirical distribution
-    [~,cdfG,distp] = palm_csort(Gvals(:),'ascend',true);
+    [~,cdfG,distp] = palm_competitive(Gvals(:),'ascend',true);
     cdfG  = unique(cdfG);
     distp = unique(distp)./numel(Gvals);
 
@@ -64,7 +64,7 @@
 else % if large G are significant (typical case)
 
     % Sort the data and compute the empirical distribution
-    [~,cdfG,distp] = palm_csort(Gvals(:),'descend',true);
+    [~,cdfG,distp] = palm_competitive(Gvals(:),'descend',true);
     cdfG  = unique(cdfG);
     distp = flipud(unique(distp))./numel(Gvals);
 

palm_lowrank.m

@@ -17,22 +17,22 @@
 
 % Common input for all cases below
 G       = varargin{1};
-[nJ,nV] = size(G);
+[nP,nV] = size(G);
 
 if nargin == 1,
 
     % Compute the basis U. This is to be done at a specific permutation.
 
     % Use a global mean, as each permutation is assumed to have the same
     % overall mean as the others. The average is the best estimate.
-    G = G - mean(G(:),1);
+    %G = G - mean(G(:),1);
 
     % Define the basis. Save some memory by working with
     % the smallest possible square of X
-    if nJ >= nV,
+    if nP >= nV,
         [U,SS,~] = svd(G'*G);
         s   = diag(SS);
-        tol = nJ * eps(max(s));
+        tol = nP * eps(max(s));
         r   = sum(s > tol);
         SSr = SS(1:r,1:r);
         U   = U(:,1:r)';
@@ -48,37 +48,66 @@
     % Outputs
     varargout{1} = diag(diag(SSr).^-.5)*U;
 
-elseif nargin == 3,
+elseif nargin == 4,
 
     if size(G,1) > 1,
 
         % With a basis known, reconstruct G, find the residuals,
         % and compute the moments.
 
         % Inputs
-        U    = varargin{2}; % basis
-        nsel = varargin{3}; % number of voxels to use
+        U     = varargin{2}; % basis
+        nsel  = varargin{3}; % number of voxels to use
+        Gmean = varargin{4}; % ensure outputs are all positive
 
         % Reconstruct the data in the new basis. Use just a subsample, so that
         % the residuals will also include this source of variability (i.e., the
         % filling of missing data, not just low-rank basis).
-        varargout{1} = zeros(size(G));
-        for p = 1:nJ,
-            idx = randperm(nV);
-            idx = idx(1:nsel);
-            varargout{1}(p,:) = G(p,idx)*pinv(U(:,idx))*U;
+        Grec = zeros(size(G));
+        if Gmean,
+            for p = 1:nP,
+                idx       = randperm(nV);
+                idx       = idx(1:nsel);
+                Grec(p,:) = (G(p,idx)-Gmean)*pinv(U(:,idx))*U + Gmean;
+                if min(Grec(p,:)) <= 0,
+                    Gfix = Grec(p,:) <= 0;
+                    Ufix = U;
+                    Ufix(:,Gfix) = -Ufix(:,Gfix);
+                    Grec(p,:) = (G(p,idx)-Gmean)*pinv(Ufix(:,idx))*Ufix + Gmean;
+                end
+            end
+        else
+            for p = 1:nP,
+                idx       = randperm(nV);
+                idx       = idx(1:nsel);
+                Grec(p,:) = G(p,idx)*pinv(U(:,idx))*U;
+            end
         end
+        varargout{1} = Grec;
 
     else
 
         % With a basis known, reconstruct a single column of G with just a few
         % known entries, and add residuals that follow a similar distribution.
 
         % Inputs
-        U    = varargin{2}; % basis
-        ysel = varargin{3}; % indices of selected tests
+        U     = varargin{2}; % basis
+        ysel  = varargin{3}; % indices of selected tests
+        Gmean = varargin{4}; % ensure outputs are all positive
 
         % Reconstruct a single permutation with lots of missing values.
-        varargout{1} = G*pinv(U(:,ysel))*U;
+        if Gmean,
+            Gmean    = Gmean;
+            Grec     = (G-Gmean)*pinv(U(:,ysel))*U + Gmean;
+            if min(Grec) <= 0,
+                Gfix = Grec <= 0;
+                Ufix = U;
+                Ufix(:,Gfix) = -Ufix(:,Gfix);
+                Grec = (G-Gmean)*pinv(Ufix(:,ysel))*Ufix + Gmean;
+            end
+        else
+            Grec = G*pinv(U(:,ysel))*U;
+        end
+        varargout{1} = Grec;
     end
 end

palm_pareto.m

@@ -71,9 +71,9 @@
     % Number of permutations & distribution CDF
     nP   = size(Gvals,1);
     if rev,
-        [~,Gvals,Gcdf] = palm_csort(Gvals,'descend',true);
+        [~,Gvals,Gcdf] = palm_competitive(Gvals,'descend',true);
     else
-        [~,Gvals,Gcdf] = palm_csort(Gvals,'ascend',true);
+        [~,Gvals,Gcdf] = palm_competitive(Gvals,'ascend',true);
     end
     Gcdf = Gcdf/nP;
 

palm_tfce.m

@@ -39,7 +39,7 @@
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
 % Choose an appropriate mask struct.
-if opts.NPC || opts.MV,
+if opts.npcmod || opts.MV,
     S = plm.maskinter;
 else
     if plm.nmasks == 1,

palm_version.txt

@@ -1 +1 @@
-Sep/2015 (version alpha58)
+Sep/2015 (version alpha59)