Refactors embedding methods and adds MCCA

Nabarb · Nabarb · commit 9cae3fec6690 · 2025-10-22T11:28:00.000+02:00
Refactors several embedding methods (PCA, GPFA, CCA, MCCA) for improved modularity and flexibility.

Adds a new MCCA embedding method, implementing multi-area canonical correlation analysis, improving multi-region analysis capabilities.

Updates argument validation for struct inputs in datareader.
diff --git a/+datareader/+is/Struct.m b/+datareader/+is/Struct.m
@@ -85,7 +85,7 @@
     fnames = fieldnames(Din);
     fnames = string(fnames);
     fields2check = ["data","time","condition","area"];
-    fieldType    = ["numeric","cell","string","string"];
+    fieldType    = ["numeric","double","string","string"];
 
     Din_ = repmat(struct(),size(Din));
 
@@ -151,11 +151,11 @@
                 'Provided time and input data dimension mismatch. Please provide,for each trial, a time struct field matching data second dimension.');
         end
     else
-        assert(iscell(opts.time),...
-            'Struct input detected. Please provide a cell array of time vectors, one per trial.');
+        % assert(iscell(opts.time),...
+        %     'Struct input detected. Please provide a cell array of time vectors, one per trial.');
 
-        assert(length(opts.time{1}) == unique(Time),...
-                'Provided time and input data dimension mismatch. Please provide a cell array of time vectors matching trial lenghts.');
+        % assert(length(opts.time{1}) == unique(Time),...
+        %         'Provided time and input data dimension mismatch. Please provide a cell array of time vectors matching trial lenghts.');
 
         dishomogeneous = false;
     end
diff --git a/+embedding/+CCA/project.m b/+embedding/+CCA/project.m
@@ -1,14 +1,16 @@
 function E = project(D,C)
 
-E = cell(size(D));
-alldata = [D{:}];
+[nTrial,nArea] = size(D);
+E = cell(nTrial,nArea);
+alldata = arrayfun(@(Aidx)cat(2,D{:,Aidx}),1:size(D,2),'UniformOutput',false);
 
-sc = (alldata - mean(alldata,2))'*C;
+E_ = cellfun(@(d,c) c * d,alldata,C,'UniformOutput',false);
 
 % For each condition, store the reduced version of each data vector
-index = 0;
-for ii = 1:length(D)
-    E{ii} = sc(index + (1:size(D(ii).data,2)),:)';
-    index = index + size(D(ii).data,2);
-end %ii
-end
+for jj = 1:nArea
+    index = 0;
+    for ii = 1:nTrial
+        E{ii,jj} = E_{jj}(index + (1:pars.TrialL),1:dims)';
+        index = index + pars.TrialL;
+    end
+end
diff --git a/+embedding/+CCA/reduce.m b/+embedding/+CCA/reduce.m
@@ -1,7 +1,7 @@
 function [E,C,Corr] = reduce(D,pars)
         % CCAREDUCE Internal function for CCA
 
-        % Agglomerate all of the conditions, and perform PCA
+        % Agglomerate all of the conditions, and perform mCCA
         E_ = cell(1,pars.nArea);
         C = cell(size(D));
         Corr = cell(size(D));
@@ -19,18 +19,19 @@
             [A,B,Corr,E_{1},E_{2}] = canoncorr(D{:});
             C = {A(:,1:dims),B(:,1:dims)};
         else
-            data = cat(1,D{:})';
+            data = arrayfun(@(Aidx)cat(2,D{:,Aidx}),1:size(D,2),'UniformOutput',false);
             d = cellfun(@(x)size(x,1),D);
-            [~,Corr,C] = embedding.CCA.mcca(data,d);
+            [~,Corr,C] = embedding.CCA.mcca(cat(1,data{:})',d);
+            E_ = cellfun(@(d,c) c * d,data,C,'UniformOutput',false);
         end
 
         % [U,V] = checkFlip(D{:},C{:},endLeg_range, interest_range);
 
         % For each condition, store the reduced version of each data vector
         E = cell(pars.nTrial,pars.nArea);
-        for ii = 1:pars.nTrial
+        for jj = 1:pars.nArea
             index = 0;
-            for jj = 1:pars.nArea
+            for ii = 1:pars.nTrial
                 E{ii,jj} = E_{jj}(index + (1:pars.TrialL),1:dims)';
                 index = index + pars.TrialL;
             end
diff --git a/+embedding/+GPFA/project.m b/+embedding/+GPFA/project.m
@@ -4,12 +4,12 @@
 alldata = [D{:}];
 projMatrix = C{:};
 
-sc = alldata'*projMatrix;
+sc = projMatrix * alldata;
 
 % For each condition, store the reduced version of each data vector
 index = 0;
 for ii = 1:length(D)
-    E{ii} = sc(index + (1:size(D{ii},2)),:)';
+    E{ii} = sc(:,index + (1:size(D{ii},2)));
     index = index + size(D{ii},2);
 end %ii
 end
diff --git a/+embedding/+GPFA/reduce.m b/+embedding/+GPFA/reduce.m
@@ -100,5 +100,5 @@
   [E{:}] = deal(AllSeq.xsm);
   [~, C{1}] = embedding.GPFA.util.orthogonalize([AllSeq.xsm], estParams.C);
   [~,lat] = pcacov(estParams.C * estParams.C');
-  VarExp{1} = cumsum(lat(1:xDim))./sum(lat);
+  VarExp{1} = cumsum(lat(xDim))./sum(lat);
 end
diff --git a/+embedding/+MCCA/loadParams.m b/+embedding/+MCCA/loadParams.m
@@ -5,25 +5,6 @@
 % pars.endLeg_range = @(t)getNormRange(t,fraction);
 % pars.interest_range = @(t)getInterestRange(t,fraction,alignment);
 % pars.ccaRefSig = [];
-pars.mcca_k = 5; 
+pars.mcca_k = 0.9; 
 
 end
-
-function T = getNormRange(t,fraction)
-    
-    Tmax = length(t);
-    T = false(size(t));
-    T(1:Tmax/fraction/2) = true;
-    T(end-Tmax/fraction/2:end) = true;
-    T = t(T);
-end
-
-function T = getInterestRange(t,fraction,alignment)
-
-    Tmax = length(t);
-    T = false(size(t));
-    Pre = alignment-Tmax/fraction/2;
-    Post = alignment+Tmax/fraction/2;
-    T(Pre:Post) = true;
-    T = t(T);
-end
diff --git a/+embedding/+MCCA/project.m b/+embedding/+MCCA/project.m
@@ -0,0 +1,16 @@
+function E = project(D,C)
+
+[nTrial,nArea] = size(D);
+E = cell(nTrial,nArea);
+alldata = arrayfun(@(Aidx)cat(2,D{:,Aidx}),1:size(D,2),'UniformOutput',false);
+
+E_ = cellfun(@(d,c) c * d,alldata,C,'UniformOutput',false);
+
+% For each condition, store the reduced version of each data vector
+for jj = 1:nArea
+    index = 0;
+    for ii = 1:nTrial
+        E{ii,jj} = E_{jj}(index + (1:pars.TrialL),1:dims)';
+        index = index + pars.TrialL;
+    end
+end
diff --git a/+embedding/+MCCA/reduce.m b/+embedding/+MCCA/reduce.m
@@ -0,0 +1,18 @@
+function [C,Corr] = reduce(D_all,pars)
+        % CCAREDUCE Internal function for CCA
+        nSets = length(D_all);
+
+        % Agglomerate all the sets, and perform mCCA
+
+        % 
+        % if length(D) < 2 && ~iscell(D)
+        %     error('Input must be a cell-array of at least two elements with data from two distinct areas.')
+        % end
+        % endLeg_range = pars.endLeg_range;
+        % interest_range = pars.interest_range;
+        
+        data = cellfun(@(eall)cat(2,eall{:}),D_all,'UniformOutput',false);
+        d = cellfun(@(x)size(x,1),data);
+        [~,Corr,C] = embedding.CCA.mcca(cat(1,data{:})',d);
+
+    end
diff --git a/+embedding/+PCA/project.m b/+embedding/+PCA/project.m
@@ -4,12 +4,12 @@
 alldata = [D{:}];
 projMatrix = C{:};
 
-sc = (alldata - mean(alldata,2))'*projMatrix;
+sc = projMatrix * (alldata - mean(alldata,2));
 
 % For each condition, store the reduced version of each data vector
 index = 0;
 for ii = 1:length(D)
-    E{ii} = sc(index + (1:size(D{ii},2)),:)';
+    E{ii} = sc(:,index + (1:size(D{ii},2)));
     index = index + size(D{ii},2);
 end %ii
 end
diff --git a/+embedding/+PCA/reduce.m b/+embedding/+PCA/reduce.m
@@ -16,11 +16,11 @@
         % For each condition, store the reduced version of each data vector
         index = 0;
         for i=1:length(D)
-            D(i).data = sc(index + (1:size(D(i).data,2)),1:dims)';
+            D(i).data = sc(index + (1:size(D(i).data,2)),:)';
             index = index + size(D(i).data,2);
         end
         [E{:}] = deal(D.data);
-        C{1} = u(:,1:dims);
-        VarExplained{1} = cumsum(lat) ./ sum(lat);  % eigenvalues
+        C{1} = u;
+        VarExplained{1} = cumsum(lat(dims)) ./ sum(lat);  % eigenvalues
 
     end
diff --git a/+embedding/GPFA.m b/+embedding/GPFA.m
@@ -1,4 +1,7 @@
-function [E,ProjMatrix,VarExplained]= GPFA(D,pars)
+function [E,ProjMatrix,VarExplained]= GPFA(D,pars,W, VarExplained)
+if nargin < 3
+    [W, VarExplained] = deal([]);
+end
 
 seqTest = pars.seqTest;
 nTrial = length(D);
@@ -21,8 +24,16 @@
 end
 [D_.T] = deal(pars.TrialL{:});
 otherArgs = [fieldnames(pars) struct2cell(pars)]';
-[E,ProjMatrix,VarExplained] = embedding.GPFA.reduce(D_(~seqTest),D_(seqTest),...
-    'xDim',pars.numPC,'verbose',false,'binWidth',pars.subsampling,...
-    otherArgs{:});
+
+if pars.projectOnly
+    E = embedding.GPFA.project(D,W);
+    ProjMatrix = W;
+else
+    [E,ProjMatrix,VarExplained] = embedding.GPFA.reduce(D_(~seqTest),D_(seqTest),...
+        'xDim',pars.numPC,'verbose',false,'binWidth',pars.subsampling,...
+        otherArgs{:});
+    ProjMatrix{1} = ProjMatrix{1}';
+end
+% VarExplained = VarExplained;
     
 end
diff --git a/+embedding/MCCA.m b/+embedding/MCCA.m
@@ -0,0 +1,24 @@
+function [E,ProjMatrix,ProjMatrixInv,Corr] = MCCA(Dall,pars)
+    nSets = numel(Dall);
+    minTrials = min(cellfun(@(x)size(x,1),Dall));
+    randIdx = arrayfun(@(idx)randi(pars(idx).nTrial,1,minTrials),1:nSets, ...
+        'UniformOutput',false);
+    Dall_ = cellfun(@(d,idx)d(idx),Dall,randIdx ...
+        ,'UniformOutput',false);
+    TrialL = arrayfun(@(p,idx)[p.TrialL{idx{1}}],pars,randIdx,'UniformOutput',false);
+
+    [ProjMatrix,Corr] = embedding.MCCA.reduce(Dall_,pars);
+    ProjMatrixInv = cellfun(@(w)pinv(w),ProjMatrix,'UniformOutput',false);
+    E_ = cellfun(@(d,c) c * cat(2,d{:}),Dall,ProjMatrix,'UniformOutput',false);
+
+
+    % For each dataset, store the coregistered version
+    E = Dall;
+    for ii = 1:nSets
+        index = 0;
+        for jj = 1:pars(ii).nTrial
+            E{ii}{jj} = E_{ii}(:,index + (1:pars(ii).TrialL{jj}),:);
+            index = index + pars(ii).TrialL{jj};
+        end
+    end
+end
diff --git a/+embedding/PCA.m b/+embedding/PCA.m
@@ -1,14 +1,26 @@
-function [E,ProjMatrix,VarExplained] = PCA(D,pars)
+function [E,ProjMatrix,ProjMatrixInv,VarExplained] = PCA(D,pars,W, VarExplained)
+
+if nargin < 3
+    [W, VarExplained] = deal([]);
+end
+
 [nT,nA] = size(D);
-E = cell(nT,nA);
-ProjMatrix = cell(1,nA);
-VarExplained = cell(1,nA);
-
-for am = 1:size(D,2)
-    D_ = repmat(struct(),size(D,1),1);
-    for id=1:size(D,1)
-        D_(id).data = full(D{id,am});
-    end
-    [E(:,am),ProjMatrix(:,am),VarExplained(:,am)] = embedding.PCA.reduce(D_,pars.numPC);
+
+D_ = repmat(struct(),size(D,1),1);
+for id=1:size(D,1)
+    D_(id).data = full(D{id});
 end
+if pars.projectOnly
+    E = embedding.PCA.project(D,W);
+    ProjMatrix = W;
+    ProjMatrixInv = {W{1}'};
+else
+    [E,ProjMatrix,VarExplained] = embedding.PCA.reduce(D_,pars.numPC);
+    E = cellfun(@(e)e(1:pars.numPC,:),E,'UniformOutput',false);
+    ProjMatrixInv = cellfun(@(w)w(:,1:pars.numPC),ProjMatrix,'UniformOutput',false);
+    ProjMatrix = cellfun(@(w)w(:,1:pars.numPC)',ProjMatrix,'UniformOutput',false);
+end
+
+
+
 end
diff --git a/+metrics/Alignment.m b/+metrics/Alignment.m
@@ -4,11 +4,12 @@
 
 pars =  metrics.pars.Alignment();
 
-theta = metrics.compute.alignment(cat(3,E{:}),pars.type,pars.nsec);
+TA = metrics.compute.alignment(cat(3,E{:}),pars.type,pars.nsec);
 
-TA = mean(...                           % mean across points
-    std(theta,[],1)...                  % for each point, compute std across trials
-    );
+% No need fot std here, already computed in the method above
+% TA = mean(...                           % mean across points
+%     std(theta,[],1)...                  % for each point, compute std across trials
+%     );                                    %
 
 end
 
diff --git a/+metrics/Tangling.m b/+metrics/Tangling.m
@@ -149,7 +149,8 @@
 %% Write everything to table
 hasDeficit = (max(succ(:,1:3),[],2)-min(succ(:,4:5),[],2) ) >.3;
 hasDeficit(6) = true;
-t = table(th(:),mixedConds(T(:))',an(:),area(:),succ(:),Deficit(:),Deficit_(:),hasDeficit(an(:)),'VariableNames',{'tangl','block','animal','area','score','deficit','deficitB','isLesioned'});
+t = table(th(:),mixedConds(T(:))',an(:),area(:),succ(:),Deficit(:),Deficit_(:),hasDeficit(an(:)), ...
+    'VariableNames',{'tangl','block','animal','area','score','deficit','deficitB','isLesioned'});
 writetable(t,'C:\Users\Fede\Documents\MATLAB\R03 figures\Results\Tangl\Tangl.GPFA.csv')
 %% New tanglind metric
 function [theta ] = newTangling(data,type,nsec)
diff --git a/@NeuralEmbedding/NeuralEmbedding.m b/@NeuralEmbedding/NeuralEmbedding.m
diff --git a/@NeuralEmbedding/computeMetrics.m b/@NeuralEmbedding/computeMetrics.m
diff --git a/@NeuralEmbedding/findEmbedding.m b/@NeuralEmbedding/findEmbedding.m
