clear all; close all; S = 0.01; BKJ = [];
 
for s = 1:400,
 
    lambda = S*(s-1);
    randn('state',0); rand('state',0); M = 2; N = 800; K = 8; e = 0.5;
    X = randn(M,N);
    Y = 0.5*randn(M,K);
    l = log2(K)*ones(1,size(Y,2));
    F = 200; BK = zeros(F,4);
 
    for i=1:F-1,
        [Y,p] = XYltoYp_MF(X,Y,l,lambda);
        Y = Y(:,find(p~=0)); 
        p = p(find(p~=0));
        p = p/sum(p); 
        % Codeword Length Update
        l = HuffLen(p);
    end; 
 
    [Y,p,D] = XYltoYpD(X,Y,l,lambda);
    p = p(find(p~=0)); p = p/sum(p);
    H = -sum(p.*log2(p));

    BKJ = [BKJ ; [lambda D H D+lambda*H size(Y,2)]]; [s lambda D H D+lambda*H size(Y,2)]
 
end;
 
plot(BKJ(:,3),BKJ(:,2),'k.'); grid on; xlabel('H (bits per vector)'); ylabel('D (MSE)');
