function[DSTsample,DSTlearn] = toy(S,xdim,ydim,T,N,iter,DSTsample)
%[DSTsample,DSTlearn] = toy(S,xdim,ydim,T,N,iter,DSTsample)
%S = number of states for discrete variables
%xdim = dimension of continuous hidden states
%ydim = dimension of data
%T = number of time steps for data (vector for multiple time series)
%N = number of time series
%iter = number of EM iterations
%DSTsample = optional argument for the DST to sample from, if not passed will use default/random DST
%
%Code to sample from a DST, learn a new one, and perform inference

if(nargin<7)
    %parent chain
	DSTsample(1).T = T; %length of each N time series (vector) 
	DSTsample(1).N = N; % number of time series
	DSTsample(1).type = 1; %type = 1 is internal node (hidden markov chain)
	DSTsample(1).S = S; %number of states
	DSTsample(1).child(1) = 2; %child number
	DSTsample(1).child(2) = 3;
	DSTsample(1).phi = normalise(rand(DSTsample(1).S,1)/20 + .5); %parameter of initial state distribution
%     DSTsample(1).Phi = mk_stochastic(rand(DSTsample(1).S,DSTsample(1).S)/20 + .5)';
    for i=1:S
        for j=1:S
            if(i==j)
        	    DSTsample(1).Phi(i,j) = .95; %parameter for state transition matrix
            else
                DSTsample(1).Phi(i,j) = .05/(S-1);
            end
        end
    end
	
	DSTsample(2).T = T;
	DSTsample(2).N = N;
	DSTsample(2).type = 1;
	DSTsample(2).S = S;
	DSTsample(2).child(1) = 4;
	DSTsample(2).parent = 1; %parent pointer
	for k=1:S
%         DSTsample(2).Phi(:,:,k) = mk_stochastic(rand(DSTsample(2).S,DSTsample(2).S)/20 + .5)';
        DSTsample(2).phi(:,k) = normalise(rand(DSTsample(2).S,1)/20 + .5,1);
        for i=1:S
            for j=1:S
                if(i==j)
            	    DSTsample(2).Phi(i,j,k) = .75;
                else
                    DSTsample(2).Phi(i,j,k) = .25/(S-1);
                end
            end
        end
	end
    
%     DSTsample(2).Phi(:,:,1) = [.95 .75; .05 .25];
%     DSTsample(2).Phi(:,:,2) = [.25 .05; .75 .95];
% 	
	DSTsample(3).T = T;
	DSTsample(3).N = N;
	DSTsample(3).type = 1;
	DSTsample(3).S = S;
	DSTsample(3).child(1) = 5;
	DSTsample(3).parent = 1;
	for k=1:S
%         DSTsample(3).Phi(:,:,k) = mk_stochastic(rand(DSTsample(3).S,DSTsample(3).S)/20 + .5)';
        DSTsample(3).phi(:,k) = normalise(rand(DSTsample(3).S,1)/20 + .5,1);
        for i=1:S
            for j=1:S
                if(i==j)
            	    DSTsample(3).Phi(i,j,k) = .75;
                else
                    DSTsample(3).Phi(i,j,k) = .25/(S-1);
                end
            end
        end
    end
    
%     DSTsample(3).Phi(:,:,1) = [.25 .05; .75 .95];
%     DSTsample(3).Phi(:,:,2) = [.95 .75; .05 .25];

	DSTsample(4).T = T;
	DSTsample(4).N = N;
	DSTsample(4).type = 3; %type 3 is a leaf node LDS
	DSTsample(4).xdim = xdim; %dimension of hidden state
	DSTsample(4).ydim = ydim; %dimension of observed emmission variable
	DSTsample(4).parent = 2;
	
% 	DSTsample(4).C = rand(ydim,xdim);
	DSTsample(4).C = eye(2); %linear operator transforming hidden state to emmission
% 	DSTsample(4).R = rand(ydim);
% 	DSTsample(4).R = DSTsample(4).R * DSTsample(4).R' + 1e-8 * eye(ydim);
%     
    DSTsample(4).R = eye(ydim)*.00001; %emmission noise covariance
	
	for i = 1:S
        
		Q = rand(xdim); %process noise covariance
		DSTsample(4).Q(:,:,i) = Q * Q' + 1e-8 * eye(xdim);
		DSTsample(4).Qinv(:,:,i) = safeinv(DSTsample(4).Q(:,:,i));
		Q1 = rand(xdim); %initial state covariance
		DSTsample(4).Q1(:,:,i) = Q1 * Q1' + 1e-8 * eye(xdim);
		DSTsample(4).Q1inv(:,:,i) = safeinv(DSTsample(4).Q1(:,:,i));
		DSTsample(4).mu(:,:,i)= rand(xdim,1); %initial state mean
	end
    
    DSTsample(4).A(:,:,1) = eye(xdim); %state transition matrix
     DSTsample(4).A(:,:,2) = [cos(pi/4) sin(pi/4); -sin(pi/4) cos(pi/4)]; 
     

	DSTsample(5).T = T;
	DSTsample(5).N = N;
	DSTsample(5).type = 3;
	DSTsample(5).xdim = xdim;
	DSTsample(5).ydim = ydim;
	DSTsample(5).parent = 3;
	
% 	DSTsample(5).C = rand(ydim,xdim);
    DSTsample(5).C = eye(2);
% 	
% 	DSTsample(5).R = rand(ydim);
% 	DSTsample(5).R = DSTsample(5).R * DSTsample(5).R' + 1e-8 * eye(ydim);
	
    DSTsample(5).R = eye(ydim)*.00001; 

	for i = 1:S
		[A1, A2, A3] = svd(rand(xdim));
        for x=1:xdim
            if( A2(x,x) > 1 )
                A2(x,x) = rand;
            elseif( A2(x,x) < -1)   
                A2(x,x) = -1 * rand;
            end
        end
        DSTsample(5).A(:,:,i) = A1 * A2 * A3';
		Q = rand(xdim);
		DSTsample(5).Q(:,:,i) = Q * Q' + 1e-8 * eye(xdim);
		DSTsample(5).Qinv(:,:,i) = safeinv(DSTsample(5).Q(:,:,i));
		Q1 = rand(xdim);
		DSTsample(5).Q1(:,:,i) = Q1 * Q1' + 1e-8 * eye(xdim);
		DSTsample(5).Q1inv(:,:,i) = safeinv(DSTsample(5).Q1(:,:,i));
		DSTsample(5).mu(:,:,i)= rand(xdim,1);
	end  
    
%     DSTsample(5).A(:,:,1) = eye(xdim) * .90;
%     DSTsample(5).A(:,:,2) = eye(xdim) * 1.10;
    
[DSTsample] = DST_sample(DSTsample);

end

DSTlearn(1).T = DSTsample(1).T;
DSTlearn(1).N = DSTsample(1).N;
DSTlearn(1).type = 1;
DSTlearn(1).S = DSTsample(1).S;
DSTlearn(1).child(1) = 2;
DSTlearn(1).child(2) = 3;

DSTlearn(2).T = DSTsample(2).T;
DSTlearn(2).N = DSTsample(2).N;
DSTlearn(2).type = 1;
DSTlearn(2).S = DSTsample(2).S;
DSTlearn(2).child(1) = 4;
DSTlearn(2).parent = 1;

DSTlearn(3).T = DSTsample(3).T;
DSTlearn(3).N = DSTsample(3).N;
DSTlearn(3).type = 1;
DSTlearn(3).S = DSTsample(3).S;
DSTlearn(3).child(1) = 5;
DSTlearn(3).parent = 1;

DSTlearn(4).T = DSTsample(4).T;
DSTlearn(4).N = DSTsample(4).N;
DSTlearn(4).type = 3;
DSTlearn(4).xdim = DSTsample(4).xdim;
DSTlearn(4).ydim = DSTsample(4).ydim;
DSTlearn(4).parent = 2;
DSTlearn(4).y = DSTsample(4).y;
for n=1:N
    
    DSTlearn(4).missing{n} = zeros(1,DSTlearn(4).T(n)); %boolean variable stating whether data is missing, if missing infers and puts inference in x
end

DSTlearn(5).T = DSTsample(5).T;
DSTlearn(5).N = DSTsample(5).N;
DSTlearn(5).type = 3;
DSTlearn(5).xdim = DSTsample(5).xdim;
DSTlearn(5).ydim = DSTsample(5).ydim;
DSTlearn(5).parent = 3;
DSTlearn(5).y = DSTsample(5).y;
for n=1:N
    
    DSTlearn(5).missing{n} = zeros(1,DSTlearn(5).T(n));
end

if(iter~=0)
 	[DSTlearn] = DST_init(DSTlearn,1,0);
	
	[DSTlearn, b] = DST_em2(DSTlearn, iter, 1,0);
end