Update MI functions

Connectivity/pairwiseMutualInformation.m

@@ -1,11 +1,15 @@
-function [I,ind,indi,indj] = pairwiseMutualInformation(X)
+function [I,ind,indi,indj] = pairwiseMutualInformation(X,normalize)
 X = X(:,:);
 n = size(X,2);
 I = zeros(n);
 ind = triu(ones(n),1);
 [indi,indj] = ind2sub([n n],find(ind(:)));
-n = length(indi);
-for k=1:n
-    I(indi(k), indj(k)) = mutualInformation(X(:,indi(k)),X(:,indj(k)));
-    I(indj(k), indi(k)) = I(indi(k), indj(k));
+nv = length(indi);
+Iv = zeros(nv);
+parfor k=1:nv
+    Iv(k) = mutualInformation(X(:,indi(k)),X(:,indj(k)),normalize);
 end
+for k=1:nv
+    I(indi(k), indj(k)) = Iv(k);
+    I(indj(k), indi(k)) = Iv(k);
+end

Connectivity/pairwiseMutualInformationROI.m

@@ -0,0 +1,34 @@
+function [I,ind,indi,indj] = pairwiseMutualInformationROI(X, hm)
+roi = hm.indices4Structure(hm.atlas.label);
+X = X(:,:);
+n = size(roi,2);
+Xc = zeros(64,size(X,1),n);
+for k=1:n
+    Xc(:,:,k) = geometricTools.resampleROI(X(:,roi(:,k))',hm.cortex.vertices, hm.cortex.faces, find(roi(:,k)), 64);
+end
+
+I = zeros(n);
+ind = triu(ones(n),1);
+[indi,indj] = ind2sub([n n],find(ind(:)));
+nv = length(indi);
+Iv = zeros(n);
+tic
+parfor k=1:nv
+%     ni = sum(roi(:,indi(k)));
+%     nj = sum(roi(:,indj(k)));
+%     if ni<nj
+%         xi = geometricTools.resampleROI(X(:,roi(:,indi(k)))',hm.cortex.vertices, hm.cortex.faces, find(roi(:,indi(k))), nj);
+%         xj = X(:,roi(:,indj(k)))';
+%     else
+%         xj = geometricTools.resampleROI(X(:,roi(:,indj(k)))',hm.cortex.vertices, hm.cortex.faces, find(roi(:,indj(k))), ni);
+%         xi = X(:,roi(:,indi(k)))';
+%     end
+    xi = Xc(:,:,indi(k));
+    xj = Xc(:,:,indj(k));
+    Iv(k) = mutualInformation(xi(:),xj(:),true);
+end
+for k=1:nv
+    I(indi(k), indj(k)) = Iv(k);
+    I(indj(k), indi(k)) = Iv(k);
+end
+toc

Connectivity/pop_srcConnectivity.m

@@ -1,4 +1,4 @@
-function EEG = pop_srcConnectivity(EEG, time, connType, normalize, neig)
+function EEG = pop_srcConnectivity(EEG, time, connType, normalize, neig, useFull)
 if nargin < 2
     time = EEG.times([1 end]);
 end
@@ -11,31 +11,60 @@
 if nargin < 5
     neig = 4;
 end
+if nargin < 6
+    useFull = false;
+end
 ind = EEG.times > time(1) & EEG.times < time(2);
 [Nroi,~, Ntrials] = size(EEG.etc.src.act(:,ind,:));
 C = zeros([Nroi Nroi, Ntrials]);
+hm = headModel.loadFromFile(EEG.etc.src.hmfile);
 
-X = permute(EEG.etc.src.act(:,ind,:),[2 1 3]);
 if strcmp(connType,'te')
+    X = permute(EEG.etc.src.act(:,ind,:),[2 1 3]);
     theta = zeros([Nroi Nroi, Ntrials]);
     for k=1:Ntrials
         disp(['Trial ' num2str(k) '/' num2str(Ntrials)]);
         [C(:,:,k), theta(:,:,k)] = conditionalTransferEntropy(X(:,:,k), -1, normalize, neig);
     end
     EEG.etc.conn = struct('measure','conditional transfer entropy','C',C,'theta',theta);
 elseif strcmp(connType,'mi')
-    [C(:,:,1),ind,indI,indJ] = pairwiseMutualInformation(X(:,:,1));
-    ind = find(ind(:)==1);
-    Cv = zeros(length(ind),Ntrials);
-    tmp = C(:,:,1);
-    Cv(:,1) = tmp(ind);
-    for k=2:Ntrials
-        disp(['Trial ' num2str(k) '/' num2str(Ntrials)]);
-        C(:,:,k) = pairwiseMutualInformation(X(:,:,k));
-        tmp = C(:,:,k);
-        Cv(:,k) = tmp(ind);
+    if useFull
+        X = permute(EEG.etc.src.actFull(EEG.etc.src.indG,ind,:),[2 1 3]);
+        if normalize
+            X = bsxfun(@rdivide, X, sqrt(sum(sum(X.^2))));
+        end
+        [C(:,:,1),ind,indI,indJ] = pairwiseMutualInformationROI(X(:,:,1), hm);
+        ind = find(ind(:)==1);
+        Cv = zeros(length(ind),Ntrials);
+        tmp = C(:,:,1);
+        Cv(:,1) = tmp(ind);
+        for k=2:Ntrials
+            disp(['Trial ' num2str(k) '/' num2str(Ntrials)]);
+            C(:,:,k) = pairwiseMutualInformationROI(X(:,:,k), hm);
+            tmp = C(:,:,k);
+            Cv(:,k) = tmp(ind);
+        end
+    else
+        X = EEG.etc.src.act(:,ind,:);
+        if normalize
+            X = bsxfun(@rdivide, X, sqrt(sum(sum(X.^2))));
+        end
+        C = pairwiseMutualInformation(X(:,:)', normalize);
+        X_surrogate = X(:,:)';
+        ntnt = size(X_surrogate,1);
+        for k=1:size(X_surrogate,2)
+            X_surrogate(:,k) = X_surrogate(randperm(ntnt,ntnt),k);
+        end
+        C_surrogate = pairwiseMutualInformation(X_surrogate, normalize);
+        X = permute(X,[2 1 3]);
+        disp(['Trial 1/' num2str(Ntrials)]);
+        [Ck(:,:,1),ind,indI,indJ] = pairwiseMutualInformation(X(:,:,1),normalize);
+        for k=2:Ntrials
+            disp(['Trial ' num2str(k) '/' num2str(Ntrials)]);
+            Ck(:,:,k) = pairwiseMutualInformation(X(:,:,k),normalize);
+        end
     end
-    EEG.etc.conn = struct('measure','mutual information','C',C,'Cv',Cv,'indI',indI,'indJ',indJ);
+    EEG.etc.conn = struct('measure','mutual information','C',C,'C_surrogate',C_surrogate,'Ck',Ck,'ind',ind,'indI',indI,'indJ',indJ);
 end
 disp('done!')