% Matlab code for simulations and experiments in IJCNN2006 paper
% "Learning to segment any random vector"
% Aapo Hyvärinen and Jukka Perkiö
% 23 March 2006
% see end of file for release notes

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% STRUCTURAL PARAMETER SETTINGS     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%parameters for data creation and estimation
subspacesize=3; %how many components in a subspace, "k" in the paper
T=1000; %number of data points
n=50; %dimension of data
sigma2=0.01^2; % noise level, must be very small

%algorithmic parameter:
maxiter=200; %how many iterations in alternating optimization estimation 

for randomdata=0:1 
%when randomdata=1, create data with completely random segmentation
%when randomdata=0, create data with Boltzmann machine with intuitive structure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CREATE DATA ACCORDING TO GENERATIVE MODEL %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~randomdata
disp('creating Boltzmann data')
     else
disp('creating random data')
end

%create model "mixing" matrix A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A=orth(randn(n,subspacesize));

%create matrix B for boltzmann machine (BM) to be used in creating S
%(random data is generated by a fake boltzmann machine with 
% parameter matrix B being all zeros )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~randomdata
%a simple beighbourhood structure
v=zeros(1,n);
v(1:floor(n/8))=1;
Bfull=toeplitz(v);
B=Bfull-diag(diag(Bfull));
else
%no structure at all
B=zeros(n);
end

%create segmentations s and observed data x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S=zeros(n,T);
X=zeros(n,T);

%create initial segmentations
for t=1:T

%discard segmentations that are just one segment
while sum(diff(S(:,t)))==0 
%(actually, this condition is too stringent because it discards 
%many other segmentations as well, but I realized this bug only 
%after final paper submission so I leave it like this

%create segmentation by gibbs sampling in BM
s=sign(randn(n,1));
for biter=1:500*(1-randomdata)
i=ceil(rand*n);
y=B(i,:)*s;
s(i)=(rand<1./(1+exp(-2*y)))*2-1;
end
S(:,t)=s;

end

%create coefficients in each segment
y_plus=randn(subspacesize,1);
y_minus=randn(subspacesize,1);

%create observed data point
X(:,t)=((S(:,t)+1)/2).*(A*y_plus)+((1-S(:,t))/2).*(A*y_minus)+randn(n,1)*sqrt(sigma2);

end % of for t



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ESTIMATE PARAMETERS                %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%INPUT DATA TO THIS PART IS IN X which contains the x(t) as its columns
%you need to set also the structural parameters, see beginning of file

%Main results are output as
%  Ahat, estimate of A
%    and 
%  Shat, which contains estimates of s(t) as its columns

disp('estimating parameters')

%initialize parameters
%%%%%%%%%%%%%%%%%%%%%%

%uniform initial value for S
Shat=ones(n,T);

%initialize A using PCA
Cx=X*X'/T;
[E, D] = eig(Cx);
[dummy,order] = sort(diag(-D));
E = E(:,order);
Ahat=E(:,1:subspacesize);

%initialize  variables
Serror=zeros(1,maxiter);
Y_plus=zeros(subspacesize,T);
Y_minus=zeros(subspacesize,T);

%START ITERATION
%%%%%%%%%%%%%%%%

for iter=1:maxiter %how many times through the data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%show progress of estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if rem(iter,10)==1
disp(maxiter-iter+1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%estimate latent variables given other parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for t=1:T %for all data points

%shortcut notation
x=X(:,t);

%estimate y_plus and y_minus for a single point, given A and s
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plusindices=find(Shat(:,t)==1);
minusindices=find(Shat(:,t)==-1);
A_plus=Ahat(plusindices,:);
A_minus=Ahat(minusindices,:);
Y_plus(:,t)=pinv(A_plus)*x(plusindices);
Y_minus(:,t)=pinv(A_minus)*x(minusindices);

if size(A_plus,1)==0, Y_plus(:,t)=zeros(subspacesize,1); end
if size(A_minus,1)==0, Y_minus(:,t)=zeros(subspacesize,1); end

%estimate s for a single point, given A and y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Shat(:,t)=sign(x.*(Ahat*(Y_plus(:,t)-Y_minus(:,t)))-.5*(Ahat*Y_plus(:,t)).^2+.5*(Ahat*Y_minus(:,t)).^2);

end %for t


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%estimate A given s and y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Basic pseudoinverse solution:

for i=1:n;

plusindices=find(Shat(i,:)==1);
minusindices=find(Shat(i,:)==-1);
Ychosen=zeros(subspacesize,T);
Ychosen(:,plusindices)=Y_plus(:,plusindices);
Ychosen(:,minusindices)=Y_minus(:,minusindices);

Ahat(i,:)=X(i,:)*pinv(Ychosen);

end % of for i

%Ahat constrained to be orthogonal
Ahat=Ahat*sqrtm(inv(Ahat'*Ahat));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end %of estimation iteration loop "for iter"


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ANALYZE AND PLOT RESULTS %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%plot matches of estimated and real segmentations

clf
hist(abs(sum(S.*Shat/n)),20)
axis('on')
axis([0,1,0,1000])

pause

%show data for the last point (because we only know its ubest)

if ~randomdata 
testpoint=1;
bar(S(:,testpoint))
axis([1,50,-2,2])
pause
bar(Shat(:,testpoint))
axis([1,50,-2,2])
pause
bar(X(:,testpoint))
axis([1,50,-2,2])
end


end %for the loop "for randomdata"

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% NOTES
%
% 1. This code does not exactly produce the simulations in the paper 
%    because I forgot to record all the random seeds needed.
%    However, there is no qualitative difference in the results.
%
% 2. The main weakness of the method that can be seen in these
%    simulations is that the noise level sigma has to be very low,
%    otherwise the results are quite bad.
%
% 3. The experiments on real data in the paper were done 
%    with essentially the same code.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
