kd_trees.m

greedy em 混和模型训练算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Kd-TREE PARTITIONING SCHEME                  %%%%%%%%%%%%%
%%%%%% Constructs m candidates for initializing     %%%%%%%%%%%%%
%%%%%% FOR DISCOVERING MOTIFS OF LENGTH W           %%%%%%%%%%%%%
%%%%%%                                                      %%%%%
%%%%%% Kostas Blekas, 19 Dec. 2001                          %%%%%
%%%%%% please contact at [e-mail] in case of problems %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [C,CP,InC] = kd_trees(X,I,C,CP,InC,b,alpha,L,W,thres)

XI = X(I,:);
ni = size(XI,1);

bgen = zeros(L,W);

for i = 1:W
    ll = zeros(L,1);
    for j = 1:L
        ll(j) = size(find(XI(:,i)==alpha(j)),1);
        if (ll(j)==0)
            ll(j) = 0.000001;
        end
    end
    bgen(:,i) = ll ./ sum(ll);
end

if (ni > 0 & ni <= b)
    [r,deg] = bestpos(XI,ni,W,thres);
end

if ( (ni > 0 & ni <= b) & (deg==1) )

   C = [C; XI(r,:)];
   
   x1 = -1 * ones(1,b);
   x1(1:ni) = I;
   InC = [InC; x1];
   
   bgen = bgen + 0.0001;
   for i=1:W
       bgen(:,i) = bgen(:,i) ./ sum(bgen(:,i));
   end
   
   CP = [CP; bgen];
else 
    k1 = repmat(max(bgen),L,1); % creates a LxW matrix with the maximum values at each column
   
    bgen1 = bgen;
    for i=1:W
        a4 = find(bgen1(:,i)<0.0001);
        bgen1(find(bgen(:,i)<0.0001),i) = k1(i);
    end
    [m1,pos]=max(sum((bgen1 - k1).^2)); % position of largest scatter
    [b1,inds] = sort(bgen(:,pos));      % sorting by frequency

    part1 = zeros(ni,1);
    part2 = zeros(ni,1);

    for i=L:-2:2
        part1(find(XI(:,pos)==alpha(inds(i))))=1;
        part2(find(XI(:,pos)==alpha(inds(i-1))))=1;
    end

    Ileft = I(find(part1==1));
    Iright = I(find(part2==1));

    % recursively call kd-tree on subpartitions
%    [C,CP,InC] = kd_recurse(X,Ileft,C,CP,InC,b,alpha,L,W,thres);
%    [C,CP,InC] = kd_recurse(X,Iright,C,CP,InC,b,alpha,L,W,thres);
    
    [C,CP,InC] = kd_recurse(X,Ileft,C,CP,InC,b,alpha,L,W,thres);
    [C,CP,InC] = kd_recurse(X,Iright,C,CP,InC,b,alpha,L,W,thres);

end