randSeed= 12;
N= 900;
N_mosaic_start=301;
N_mosaic_end=600;
tt_ratio=2;
RV= 0.2; % Factor of rate variation

% Branch lengths
w0=0.15;	
w=[w0,w0,w0,w0,2*w0];

% ---------------------------------------
%  Random number generator seed
% ---------------------------------------
rand('state',randSeed);

% ---------------------------------------
%          	 Hidden nodes 
% ---------------------------------------
for i=1:N
	hidden(1,i)=SampleFromMultinomial([0.25 0.25 0.25 0.25]);
end

for i=1:N_mosaic_start-1
	hidden(2,i)=Mutation(hidden(1,i),w(5),tt_ratio);
end
for i=N_mosaic_end+1:N
	hidden(2,i)=Mutation(hidden(1,i),w(5),tt_ratio);
end

% ---------------------------------------
%     Leave nodes 
% ---------------------------------------
for i=1:N_mosaic_start-1
	species(1,i)=Mutation(hidden(1,i),w(1),tt_ratio);
	species(2,i)=Mutation(hidden(1,i),w(2),tt_ratio);
	species(3,i)=Mutation(hidden(2,i),w(3),tt_ratio);
	species(4,i)=Mutation(hidden(2,i),w(4),tt_ratio);
end
for i=N_mosaic_end+1:N
	species(1,i)=Mutation(hidden(1,i),w(1),tt_ratio);
	species(2,i)=Mutation(hidden(1,i),w(2),tt_ratio);
	species(3,i)=Mutation(hidden(2,i),w(3),tt_ratio);
	species(4,i)=Mutation(hidden(2,i),w(4),tt_ratio);
end


% ---------------------------------------
%     Region subject to rate heterogeneity
% ---------------------------------------

for i=N_mosaic_start:N_mosaic_end
	hidden(2,i)=Mutation(hidden(1,i),RV*w(5),tt_ratio);
end

for i=1:N
	species(1,i)=Mutation(hidden(1,i),RV*w(1),tt_ratio);
	species(2,i)=Mutation(hidden(1,i),RV*w(2),tt_ratio);
	species(3,i)=Mutation(hidden(2,i),RV*w(3),tt_ratio);
	species(4,i)=Mutation(hidden(2,i),RV*w(4),tt_ratio);
end

Data(1:4,1:N,1)= species;

save Data.mat Data;


dummy=input('The data have been saved in Data.mat. Press any key to finish.');