%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%    demonstration of EM using binary images               %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clf;
 
% load images

image1=[
  0 0 0 0 0 0 0 0 0
  0 0 0 0 0 0 0 0 0
  0 0 1 1 1 1 1 0 0
  0 0 1 0 0 0 1 0 0 
  0 0 1 0 0 0 1 0 0 
  0 0 1 0 0 0 1 0 0 
  0 0 1 1 1 1 1 0 0 
  0 0 0 0 0 0 0 0 0 
  0 0 0 0 0 0 0 0 0];

image2=[
  0 0 0 0 0 0 0 0 0 
  0 1 0 0 0 0 0 1 0 
  0 0 1 0 0 0 1 0 0 
  0 0 0 1 0 1 0 0 0 
  0 0 0 0 1 0 0 0 0 
  0 0 0 1 0 1 0 0 0 
  0 0 1 0 0 0 1 0 0 
  0 1 0 0 0 0 0 1 0 
  0 0 0 0 0 0 0 0 0];

image3=[
  0 0 0 0 0 0 0 0 0 
  0 1 1 1 1 1 1 1 0 
  0 1 1 1 1 1 1 1 0 
  0 1 1 1 1 1 1 1 0 
  0 1 1 1 1 1 1 1 0 
  0 1 1 1 1 1 1 1 0 
  0 1 1 1 1 1 1 1 0 
  0 1 1 1 1 1 1 1 0 
  0 0 0 0 0 0 0 0 0];

image4=[
  1 1 1 1 1 1 1 1 1 
  1 1 1 1 1 1 1 0 1 
  1 1 1 1 1 1 0 0 1 
  1 1 1 1 1 0 1 0 1 
  1 1 1 1 0 1 1 0 1 
  1 1 1 0 1 1 1 0 1 
  1 1 0 1 1 1 1 0 1 
  1 0 0 0 0 0 0 0 1 
  1 1 1 1 1 1 1 1 1];

image5=[
  0 0 0 0 0 0 0 0 0 
  0 0 0 0 1 0 0 0 0 
  0 0 0 0 1 0 0 0 0 
  0 0 0 0 1 0 0 0 0 
  0 1 1 1 1 1 1 1 0 
  0 0 0 0 1 0 0 0 0 
  0 0 0 0 1 0 0 0 0 
  0 0 0 0 1 0 0 0 0 
  0 0 0 0 0 0 0 0 0];

% given parameters

[lx,ly]=size(image1);
num = lx*ly;

% defining parameters

P   = 0.2;
pm  = .3;

M   = 5;

ITER= 20;
T   = 40;
alpha=.9;

rand('seed',65875);

% setting up images

world(1,:)=reshape(image1,1,num);
world(2,:)=reshape(image2,1,num);
world(3,:)=reshape(image3,1,num);
world(4,:)=reshape(image4,1,num);
world(5,:)=reshape(image5,1,num);

cell=round(rand(M,num));

% show cell initial values

for m=1:M,
  subplot(3,M,m);
  colormap(gray)
  image(256*reshape(cell(m,:),lx,ly));
  axis equal; axis off;
end

% show existing outer world model

for m=1:M,
  subplot(3,M,2*M+m);
  colormap(gray)
  image(256*reshape(world(m,:),lx,ly));
  axis equal; axis off;
end

% main roop

for iter=1:ITER,
  
  for t=1:T,
    tmp=randperm(M);
    index(t)=tmp(1);
  end

  data=xor(world(index,:),(round(rand(T,num)+P-1/2)));
  
  % show data
  
  for t=1:T,
    subplot(3,M,ceil(3*M/2));
    colormap(gray)
    image(256*reshape(data(t,:),lx,ly));
    axis equal; axis off;
    drawnow
  end
  
  for t=1:T,
    tmp=ones(M,1)*data(t,:);
    diff=xor(tmp,cell);
    hamming(:,t)=diff*ones(num,1);
    p(:,t)=(pm.^(hamming(:,t))).*((1-pm).^(num-hamming(:,t)));
  end

  psum=sum(p);
  pcon=p./(ones(M,1)*psum);
  
  average=zeros(M,num);
  for t=1:T,
    average=average+pcon(:,t)*data(t,:);
  end

  average=average./(pcon*ones(T,num));
  
  cell=round((1-alpha)*cell+alpha*average);
   
  % show cell values

  for m=1:M,
    subplot(3,M,m);
    colormap(gray)
    image(256*reshape(cell(m,:),lx,ly));
    axis equal; axis off;
  end
  drawnow

  fprintf('hit return\n');
  pause;
  
  for t=1:T,
    tmp=ones(M,1)*data(t,:);
    diff=xor(tmp,cell);
    hamming(:,t)=diff*ones(num,1);
  end
  pm=ones(1,M)*(hamming.*pcon)*ones(T,1)/(T*num);
end