add

AlginPhi.m

@@ -0,0 +1,21 @@
+function [perm_phi,perm_idx] = AlginPhi(basephi, phi)
+
+    [D,K] = size(basephi.betaMu);
+    perm_phi = struct('alpha',zeros(1,K),'v',D*ones(1,K),'invWBetaMuMu',zeros(D,D,K),...
+                'betaMu',zeros(D,K),'beta',zeros(1,K));
+
+    basemu = zeros(D,K);    
+    mu = basemu;
+    for i = 1:K
+       basemu(:,i) = basephi.betaMu(:,i)/basephi.beta(i);   
+       mu(:,i) = phi.betaMu(:,i)/phi.beta(i);   
+    end 
+
+     perm_idx = cluster_algin( mu', basemu' );
+
+     perm_phi.alpha = phi.alpha(:,perm_idx);
+     perm_phi.v = phi.v(:,perm_idx);
+     perm_phi.invWBetaMuMu = phi.invWBetaMuMu(:,:,perm_idx);
+     perm_phi.betaMu = phi.betaMu(:,perm_idx);
+     perm_phi.beta = phi.beta(:,perm_idx);
+end

CreateNetworksFunc.m

@@ -0,0 +1,41 @@
+function [ Network ] = CreateNetworksFunc(Conf)
+%CREATENETWORKS Summary of this function goes here
+%   Detailed explanation goes here
+
+    num = Conf.NodeNumber;
+    square = Conf.Square;
+    maxDist = Conf.CommDist;
+
+    loc = square*rand(num,2) - square/2;        
+
+    Dists = Euclid_Dist(loc(:,1),loc(:,2));
+
+    % without self-loop
+    Dists = Dists + 10*maxDist*eye(num);
+
+    Neighbors = cell(num,1);
+    maxDegree = 0;
+    edges = 0;
+    for i=1:num
+        Neighbors{i} = find(Dists(i,:)<=maxDist);
+        if length(Neighbors{i}) > maxDegree
+            maxDegree = length(Neighbors{i});
+        end
+        edges = edges + length(Neighbors{i});
+    end
+
+    Nodes.loc = loc;
+    Nodes.neighbors = Neighbors;
+
+    Network.maxDegree = maxDegree;
+    Network.edges = edges/2; %% undirected graph
+    Network.Conf = Conf;
+    Network.Nodes = Nodes;
+end
+
+function dist = Euclid_Dist(X,Y)
+    len = length(X);
+    xx = repmat(X,1,len);
+    yy = repmat(Y,1,len);    
+    dist = sqrt((xx-xx').^2+(yy-yy').^2);
+end

Diff_NaturalParameter.m

@@ -0,0 +1,23 @@
+function phi = Diff_NaturalParameter(idx, newphi, oldphiGroup, Neighbors, w)          
+%        
+    [D,K] = size(newphi.betaMu);
+    phi = struct('alpha',zeros(1,K),'v',D*zeros(1,K),'invWBetaMuMu',zeros(D,D,K),...
+         'betaMu',zeros(D,K),'beta',zeros(1,K));      
+
+    phi.alpha = w(idx)*newphi.alpha;
+    phi.v = w(idx)*newphi.v;
+    phi.beta = w(idx)*newphi.beta;    
+    phi.betaMu = w(idx)*newphi.betaMu;    
+    phi.invWBetaMuMu = w(idx)*newphi.invWBetaMuMu;    
+
+    for i=1:length(Neighbors)
+        nei = Neighbors(i);
+        phi.alpha = phi.alpha + w(nei)*oldphiGroup(nei).alpha;
+        phi.v =  phi.v + w(nei)*oldphiGroup(nei).v;
+        phi.beta = phi.beta + w(nei)*oldphiGroup(nei).beta;
+        phi.betaMu = phi.betaMu + w(nei)*oldphiGroup(nei).betaMu;
+        phi.invWBetaMuMu = phi.invWBetaMuMu + w(nei)*oldphiGroup(nei).invWBetaMuMu;
+    end       
+end
+
+

DrawNetworks.m

@@ -0,0 +1,25 @@
+function fig = DrawNetworks( Network )
+%DRAWNETWORKS Summary of this function goes here
+%   Detailed explanation goes here
+
+num = Network.Conf.NodeNumber;
+%maxDegree = Network.maxDegree;
+loc = Network.Nodes.loc;
+square = Network.Conf.Square;
+Neighbors = Network.Nodes.neighbors;
+
+fig = figure;
+plot(loc(:,1),loc(:,2),'ro','MarkerSize',8,'LineWidth',2);
+axis([-square/2,square/2,-square/2,square/2]);   
+for i=1:num
+    for k = 1:length(Neighbors{i})
+        j = Neighbors{i}(k);
+        %     c = num2str(Dists(i,j),'%.2f');
+        %     text((loc(i,1) + loc(j,1))/2,(loc(i,2) + loc(j,2))/2,c,'Fontsize',10);
+        %     hold on;
+        line([loc(i,1),loc(j,1)],[loc(i,2),loc(j,2)],'LineWidth',0.8,'Color','b');
+    end
+end
+set(gcf, 'Color', 'w');
+end
+

GradientDescent.m

@@ -0,0 +1,14 @@
+function phi = GradientDescent(oldphi,phi_star,eta)
+%
+    phi = oldphi;
+
+    phi.alpha = (1-eta)*phi.alpha + eta*phi_star.alpha;
+    phi.v =  (1-eta)*phi.v + eta*phi_star.v;
+    phi.beta = (1-eta)*phi.beta + eta*phi_star.beta;    
+    phi.betaMu =(1-eta)*phi.betaMu + eta*phi_star.betaMu;
+    phi.invWBetaMuMu =(1-eta)*phi.invWBetaMuMu + eta*phi_star.invWBetaMuMu;
+
+
+end
+
+

InitVB.m

@@ -0,0 +1,17 @@
+
+function R0 = InitVB(data,K)
+%
+    [N,~] = size(data);                
+    C = NaN;
+%     while(isnan(C))
+%         %[IDX,C] = kmeans(data,K,'start','uniform','emptyaction','drop');  %,
+%         [IDX, C] = litekmeans(data, K,'Replicates',20);
+%     end
+
+   [~,IDX] = max(rand(N,K),[],2);
+
+    R0 = zeros(N,K);
+    for i = 1:K
+        R0(:,i) = (IDX == i);
+    end        
+end

Metropolis_Weight.m

@@ -0,0 +1,14 @@
+function w = Metropolis_Weight(NodeNum,Neighbors)
+    % update node i
+    w = zeros(NodeNum,NodeNum);
+
+    for i=1:NodeNum
+        degree = length(Neighbors{i});                
+        for j=1:length(Neighbors{i})
+            nei = Neighbors{i}(j);
+            degree2 = length(Neighbors{nei});                
+            w(nei,i) = 1/(1+ max([degree,degree2]));
+        end        
+        w(i,i) = 1 - sum(w(:,i));
+    end
+end

MutualInfo.m

@@ -0,0 +1,57 @@
+function MIhat = MutualInfo(L1,L2)
+%   mutual information
+%
+%   version 2.0 --May/2007
+%   version 1.0 --November/2003
+%
+%   Written by Deng Cai (dengcai AT gmail.com)
+%===========    
+L1 = L1(:);
+L2 = L2(:);
+if size(L1) ~= size(L2)
+    error('size(L1) must == size(L2)');
+end
+
+Label = unique(L1);
+nClass = length(Label);
+
+Label2 = unique(L2);
+nClass2 = length(Label2);
+if nClass2 < nClass
+     % smooth
+     L1 = [L1; Label];
+     L2 = [L2; Label];
+elseif nClass2 > nClass
+     % smooth
+     L1 = [L1; Label2];
+     L2 = [L2; Label2];
+end
+
+
+G = zeros(nClass);
+for i=1:nClass
+    for j=1:nClass
+        G(i,j) = sum(L1 == Label(i) & L2 == Label(j));
+    end
+end
+sumG = sum(G(:));
+
+P1 = sum(G,2);  P1 = P1/sumG;
+P2 = sum(G,1);  P2 = P2/sumG;
+if sum(P1==0) > 0 || sum(P2==0) > 0
+    % smooth
+    error('Smooth fail!');
+else
+    H1 = sum(-P1.*log2(P1));
+    H2 = sum(-P2.*log2(P2));
+    P12 = G/sumG;
+    PPP = P12./repmat(P2,nClass,1)./repmat(P1,1,nClass);
+    PPP(abs(PPP) < 1e-12) = 1e-12; %
+    MI = sum(P12(:) .* log2(PPP(:)));
+    MIhat = MI / max(H1,H2);
+    %%%%%%%%%%%%%   why complex ?       %%%%%%%%
+    MIhat = real(MIhat);
+end
+
+
+

NormalizeFea.m

@@ -0,0 +1,60 @@
+function fea = NormalizeFea(fea,row)
+% if row == 1, normalize each row of fea to have unit norm;
+% if row == 0, normalize each column of fea to have unit norm;
+%
+%   version 3.0 --Jan/2012 
+%   version 2.0 --Jan/2012 
+%   version 1.0 --Oct/2003 
+%
+%   Written by Deng Cai (dengcai AT gmail.com)
+%
+
+if ~exist('row','var')
+    row = 1;
+end
+
+if row
+    nSmp = size(fea,1);
+    feaNorm = max(1e-14,full(sum(fea.^2,2)));
+    fea = spdiags(feaNorm.^-.5,0,nSmp,nSmp)*fea;
+else
+    nSmp = size(fea,2);
+    feaNorm = max(1e-14,full(sum(fea.^2,1))');
+    fea = fea*spdiags(feaNorm.^-.5,0,nSmp,nSmp);
+end
+
+return;
+
+
+
+
+
+
+
+if row
+    [nSmp, mFea] = size(fea);
+    if issparse(fea)
+        fea2 = fea';
+        feaNorm = mynorm(fea2,1);
+        for i = 1:nSmp
+            fea2(:,i) = fea2(:,i) ./ max(1e-10,feaNorm(i));
+        end
+        fea = fea2';
+    else
+        feaNorm = sum(fea.^2,2).^.5;
+        fea = fea./feaNorm(:,ones(1,mFea));
+    end
+else
+    [mFea, nSmp] = size(fea);
+    if issparse(fea)
+        feaNorm = mynorm(fea,1);
+        for i = 1:nSmp
+            fea(:,i) = fea(:,i) ./ max(1e-10,feaNorm(i));
+        end
+    else
+        feaNorm = sum(fea.^2,1).^.5;
+        fea = fea./feaNorm(ones(1,mFea),:);
+    end
+end
+
+

Uniform_Weight.m

@@ -0,0 +1,14 @@
+
+function w = Uniform_Weight(NodeNum,Neighbors)
+% update node i            
+    w = zeros(NodeNum,NodeNum);
+
+    for i =1:NodeNum
+        degree = length(Neighbors{i})+1;      
+        for j=1:length(Neighbors{i})
+            nei = Neighbors{i}(j);
+            w(nei,i) = 1/degree;
+        end
+        w(i,i) = 1/degree;        
+    end
+end

UpdateDualVariables.m

@@ -0,0 +1,9 @@
+function Dual = UpdateDualVariables(Dual,Phi1,Phi2,rhos)
+
+    Dual.alpha = Dual.alpha + rhos(1)/2*(Phi1.alpha-Phi2.alpha);    
+    Dual.v = Dual.v + rhos(2)/2*(Phi1.v-Phi2.v);
+    Dual.invWBetaMuMu = Dual.invWBetaMuMu + rhos(3)/2*(Phi1.invWBetaMuMu-Phi2.invWBetaMuMu);
+    Dual.betaMu = Dual.betaMu + rhos(4)/2*(Phi1.betaMu-Phi2.betaMu);
+    Dual.beta = Dual.beta + rhos(5)/2*(Phi1.beta-Phi2.beta);
+
+end

VBUpdate.m

@@ -0,0 +1,48 @@
+function [Hypers,History,R,Label] = VBUpdate(data,K,GroundTruth)
+
+    % data: N x D 
+    [~,D] = size(data);
+    Prior = struct('alpha0',1e-5,'beta0',1e-5,'mu0',1e-5*ones(D,1),'v0',D+1,'invW0',1e-5*eye(D,D));    
+
+    L0 = -inf;        
+    eps =1e-12;
+
+    maxIters = 300;
+
+    % history 
+    History.kl_error = zeros(1,maxIters);
+    History.old_Params = repmat(struct('mix',zeros(1,K),'Mu',zeros(D,K),'Sigma',zeros(D,D,K)),maxIters,1);
+
+
+    R = InitVB(data,K);
+
+    for t=1:maxIters           
+        Hypers = VBM_step(data,R,Prior);
+        [R,logR] = VBE_step(data,Hypers);
+        [Hypers,R] = AlignVBResults(Hypers,R);  
+        L = VBbound(data, Hypers,R,logR,Prior);
+
+        History.Lseq(t) = L;
+        History.kl_error(t) = kl_mixgaus(Hypers, GroundTruth);            
+        History.old_Params(t) = Hyper2Params(Hypers);        
+
+%         fprintf('%e %e \n',abs(L-L0),eps*abs(L));
+%         if abs(L-L0) < eps*abs(L)
+%             break;
+%         end          
+        L0 = L;
+        mixx = History.old_Params(t).mix;
+        fprintf('%d: %f %f %f %f\n',t,mixx(1),mixx(2),mixx(3),History.kl_error(t));
+    end
+
+    [Hypers,R] = AlignVBResults(Hypers,R);
+    [~,Label] = max(R,[],2); 
+
+
+   disp(['Total iteratons:',num2str(t)]);    
+end
+
+
+
+
+

VBfunctions/AlignVBResults.m

@@ -0,0 +1,18 @@
+function [newHypers,newR] = AlignVBResults(oldHypers,oldR,base_align)
+    %reorder the clusters
+
+    [D,K] = size(oldHypers.Mu);
+
+    newHypers = struct('invW',zeros(D,D,K),'alpha',zeros(1,K),'v',D*ones(1,K),...
+                'beta',zeros(1,K),'Mu',zeros(D,K));  
+    newR = zeros(size(oldR));        
+
+    [idxmap,cost] = cluster_algin(oldHypers.Mu', base_align');      
+
+    newHypers.alpha = oldHypers.alpha(idxmap);
+    newHypers.beta = oldHypers.beta(idxmap);  
+    newHypers.Mu = oldHypers.Mu(:,idxmap);
+    newHypers.v = oldHypers.v(idxmap);         
+    newHypers.invW = oldHypers.invW(:,:,idxmap);
+    newR = oldR(:,idxmap);
+end

VBfunctions/GaussPDF.m

@@ -0,0 +1,12 @@
+function [ prob ] = GaussPDF(Data,M,Sigma)
+%  Calculate Gaussian Probability Distribution Function
+%   Data  dim x N
+%   M     dim x 1
+%   Sigma dim x dim
+
+    [dim,N] = size(Data);
+    Data = Data'-repmat(M',N,1);
+    prob = sum((Data/Sigma).*Data,2); % Data*inv(Sigma)
+    prob = exp(-0.5*prob)/sqrt((2*pi)^dim*(abs(det(Sigma))+realmin));
+
+end