%% Markov chain simulator
% produces gif, avi(no sound) and wav(sound)
% if you want to combine animation and sound, you can use iMovie (mac), Avs
% video editor (win), ?? (linux) - if you find a better or free-er
% solution, please let me know (especially in linux)
%
% avi-writing seems to not work in Mac right now, but it worked fine in
% Windows (replace line 168 by line 167).
%
% Feel free to reuse/rewrite/redistribute this file.
%
% Alexey Morozov morozova }at} umich dot edu 19 Nov 2012
%
% References
% help from http://www.haverford.edu/econ/econ365/Note%20on%20Markov%20Chains.pdf
% stochastic matrices in http://iopscience.iop.org/0305-4470/36/12/333/
% rand matr http://sites.poli.usp.br/pmr/ltd/People/jside/IdeCozman_sbia02.pdf
% google matrix http://arxiv.org/pdf/1106.6215v2.pdf
% tone frequencies from http://www.phy.mtu.edu/~suits/notefreqs.html
% embed mp3 http://www.w3schools.com/html/html_sounds.asp


clear
clc
format compact
% format long
cd ~/Desktop/hw

n=8; %number of states

f=[261.63;293.66;329.63;349.23;392.00;440;493.88;523.25;587.33]; %C4 to C5 from %http://www.phy.mtu.edu/~suits/notefreqs.html
if n>9
    f=linspace(f(1),f(end),n);
end
nsec=0.25;%how many seconds to play each tone
stop=100; %how many steps to take
Fs=44100; %sampling frequency


x=1:n; %state space

%% simple n=3
% P=[.1 .9 0;    0 .3 .7;    .2 0 .8]; %state transition matrix

%% random
% P=rand(n);

%% 4-2-2 setup

P=zeros(n);

p1=ones(4); %first block, call it p1
p2=ones(2); %second
p3=ones(2); %third

P(1:4,1:4)=p1;
P(5:6,5:6)=p2;
P(7:8,7:8)=p3;
% P(1,n)=.5; P(n,1)=.5;
% P(6,7)=.5; P(7,6)=.5; %everybody is connected
i=.5;
P(6,2)=i;
P(8,5)=i; % before this line, there was no way to get out of p3 or p2
P(2,6)=i; % before this line, there was no way to get out of       p2
P(4,8)=i; % now, can get out of p3, p2,

%% 3-3-3

% P=zeros(n);
% p1=ones(3);
% p2=ones(3);
% p3=ones(3);
%
% P(1:3,1:3)=p1;
% P(4:6,4:6)=p2;
% P(7:9,7:9)=p3;
% % P(1,n)=.5; P(n,1)=.5;
% % P(6,7)=.5; P(7,6)=.5; %everybody is connected
% i=.5;
% P(5,2)=i;
% P(8,5)=i; %no way to get out of p3 or p2
% P(2,6)=i; %no way to get out of     p2
% P(4,8)=i; %all connected


%% chords
%C major triad, as in http://en.wikipedia.org/wiki/Chord_(music)#Triads

% P=zeros(n);
% [x,y]=meshgrid([1 3 5]);
% i=sub2ind([n n], x(:), y(:));
% P(i)=1;
%
%
% [x,y]=meshgrid([2 4 6]);
% i=sub2ind([n n], x(:), y(:));
% P(i)=1;
%
% i=.5;
% P(6,5)=i;
% P(5,6)=i;

%% all variants need this to be ran

P = bsxfun(@rdivide,P,sum(P,2)); %normalize all the rows so that they sum to 1
pi0=rand(1,n); pi0 = pi0/sum(pi0); %initial prob distr

r=rand(stop,1); %generate the random numbers which are later used to decide which state to go to next. If you don't rerun this cell every time, you'll get the same simulation

% colormap
%colormap tester
% h=jet; scatter3(h(:,1),h(:,2),h(:,3),100,h)
% h=hot; h=h(:,[2 1 3]); scatter3(h(:,1),h(:,2),h(:,3),100,h)
i=32; h=[zeros(i,1) (1:i)'/i zeros(i,1); (1:i)'/i ones(i,1) (1:i)'/i]; %custom colormap, to see what it is in RGB space, uncomment the next line - black to green to white
% scatter3(h(:,1),h(:,2),h(:,3),100,h)

%% no random number generation below this line - so if you just run the cells below this line, you'll keep getting the same trajectory

%% simulate, write gif, avi, wav, png (just P values), another png is created in the next section (P values and evolution history)

fname='mc'; %this will be the first part of all files (extensions appended later)

figure(1); clf
set(gcf,'color',[1 1 1]) %I prefer white background in my GIF's, but you can change it

imagesc(P);
ylabel('Current state')
xlabel('Next possible state')

cb=colorbar;
ylabel(cb,'Transition probability')

% colormap(h)
colormap(flipud(h)) %matrix, alhough maybe google paper was not flipud
set(gca,'clim',[0 1])

[x,y]=meshgrid(1:n);
text(x(:)-.25,y(:),num2str(round(P(:)*100)/100'),'fontsize',8)

set(gcf,'PaperPositionMode','auto')
print(gcf,'-dpng','-r0', [fname '0.png'])


hold on

c=jet(stop);
notes=[];
s=nan(1,stop);
for i=1:stop
    if i==1;
        
        
        s(1)=find(r(1)<=cumsum(pi0),1); %pick the initial state using initial distribution pi0
        
        
        x=line([s(i) s(i)],[s(i) s(i)],'color','b','linewidth',2); %line showing where the state will go next in x direction
        y=line([s(i) s(i)],[s(i) s(i)],'color','b','linewidth',2); %line showing where the state will go next in y direction
        d=plot([s(i) s(i)],[s(i) s(i)],'ro','markerfacecolor','r','markersize',10); %the current state dot
        % set([x y d],'EraseMode','none');
        t=text(n,n+1.3,['TimeStep = ' num2str(i) '/' num2str(stop) ]);
        fr=getframe(gcf);
        [im1,ma] = frame2im(fr);
        [im2,ma]=rgb2ind(im1,128);
        imwrite(im2,ma,[fname '.gif'],'LoopCount',0,'DelayTime',1);%write gif
        %         aviobj = avifile('mc.avi','colormap',ma,'compression','MSVC','quality',75,'fps',1/nsec); % Windows
        aviobj = avifile([fname '.avi'],'colormap',ma,'compression','none','fps',1/nsec); %Mac, write avi
        aviobj = addframe(aviobj,im2); %write avi
    
    else
        
        
        s(i)=find(r(i)<=cumsum( P(s(i-1),:) ),1); %pick new state at the current step using Markov Chain transition probability matrix P
        
        
        line([s(i-1) s(i) s(i)] +(i-stop/2)/(5*stop),[s(i-1) s(i-1) s(i)]+(i-stop/2)/(5*stop),'color',c(i,:),'linewidth',2) %new transition
        delete(d) %erase old red dot
        d=plot(s(i)+(i-stop/2)/(5*stop),s(i)+(i-stop/2)/(5*stop),'ro','markerfacecolor','r','markersize',10); %the current state dot
        delete(t) %erase old text
        t=text(n,n+1.3,['TimeStep = ' num2str(i) '/' num2str(stop) ]); %new text
        fr=getframe(gcf);
        im1 = frame2im(fr);
        im2=rgb2ind(im1,ma);
        imwrite(im2,ma,[fname '.gif'],'WriteMode','append'); %write gif
        aviobj = addframe(aviobj,im2); %write avi
        
    end
    
    Ws=2*pi*f(s(i))/Fs; %sampled cyclic note frequency
    
    %fixed time (nsec) of seconds (good timing for video, but produces jumps in sound)
    % note=sin( Ws*(1:(Fs*nsec) ) );
    
    %fixed number of periods
    %     np=400;
    %     note=sin( Ws*(1:(np*2*pi/Ws)) ); %add np periods (around a second)of the sinusoid at that frequency and make sure we take integer number of full periods, so there is no jumps
    
    %fixed time and number of periods (sounds smooth, but possibly introduces about 0.001 second lag or lead per note)
        ks=Fs*nsec; %how many steps nsec will take for this particular frequency
        Ts=2*pi/Ws; %period for the current note (steps)
        npr=round(ks/Ts); %what is the closest integer number of periods to nsec?
        %disp( num2str(npr/f(s(i))-nsec,'%1.4f') ) %display the synchronization error in seconds
        note=sin( Ws*(1:(npr*Ts)) ); %take the closest number of periods of current note
    
    %fixed time (np/400) of a second , multiple notes
    %     np=400;
    %     if s(i)==1;
    %     note=sin( 2*pi*f(1)/Fs*(1:(Fs*np/400) ) )+...
    %         sin( 2*pi*f(3)/Fs*(1:(Fs*np/400) ) )+...
    %         sin( 2*pi*f(5)/Fs*(1:(Fs*np/400) ) );
    %     elseif s(i)==2
    %     note=sin( 2*pi*f(2)/Fs*(1:(Fs*np/400) ) )+...
    %         sin( 2*pi*f(4)/Fs*(1:(Fs*np/400) ) )+...
    %         sin( 2*pi*f(7)/Fs*(1:(Fs*np/400) ) );
    %     else
    %     note=sin( 2*pi*f(3)/Fs*(1:(Fs*np/400) ) )+...
    %         sin( 2*pi*f(5)/Fs*(1:(Fs*np/400) ) )+...
    %         sin( 2*pi*f(8)/Fs*(1:(Fs*np/400) ) );
    %     end
    
    
    %     wavplay(note,Fs) % uncomment if you want to hear what got added right when it got added inside the loop
    
    notes=[notes note]; %store the tones
    
end

hold off

aviobj = close(aviobj); %close avi

% stairs(note(end-1000:end)) %plot the sound
% wavplay(notes,Fs) %play the whole sound
wavwrite(notes,Fs,'mc.wav') %write sound


%% history color plot in png (evolution of MC)
figure(2); clf
set(gcf,'color',[1 1 1])
imagesc(P)
ylabel('Current note')
xlabel('Next possible note')

ns={'C4','D4','E4','F4','G4','A4','B4','C5','D5'};
set(gca,'YTickLabel',ns(1:n));
set(gca,'XTickLabel',ns(1:n));
% set(gca,'YTickLabel',[]);
% yt=get(gca,'YTick');
% text(0*yt,yt,{'C_4','D_4','E_4','F_4','G_4','A_4','B_4','C_5'});
%
% set(gca,'XTickLabel',[]);
% xt=get(gca,'XTick');
% text(xt,0*xt+n+1,{'C_4','D_4','E_4','F_4','G_4','A_4','B_4','C_5'});

cb=colorbar;
ylabel(cb,'Transition probability')

% colormap(h)
colormap(flipud(h)) %matrix, alhough maybe google paper was not flipud
set(gca,'clim',[0 1])
[x,y]=meshgrid(1:n);
text(x(:)-.25,y(:),num2str(round(P(:)*100)/100'),'fontsize',8)

ss=[s;s]; %2 rows on top of each other [1 2 3...; 1 2 3...]
ss=ss(:); %1 column with [1 1 2 2 3 3 ...]
x=ss(2:end);
y=ss(1:(end-1) );

hold on
c=flipud(gray(stop));
c=jet(stop);
for i=1:stop-1
    %   line(x(2*i-1+(0:2)),y(2*i-1+(0:2)),'color',c(i,:),'linewidth',3*(stop-i+1)/stop)
    line(x(2*i-1+(0:2))+(i-stop/2)/(5*stop),y(2*i-1+(0:2))+(i-stop/2)/(5*stop),'color',c(i,:),'linewidth',2)
end
hold off
set(gcf,'PaperPositionMode','auto')
print(gcf,'-dpng','-r0','mc.png')

%% limiting distr (time spent) (can we go backward (generate P based on pii) ?)
% [V,L]=eig(P);
% Q=L^100; %only the evalue=1 persists
% W=V^(-1);
% W1=abs(W(1,:));
% W1=W1/sum(W1) %limiting distr
% 
% %another way:
% % pii*P=pii %limiting eqn
% % pii*(P-I)=0 %limiting eqn
% % (P-I)'*pii'=0 %limiting eqn
% % inv((P-eye(n))')*zeros(n,1) %gets 0
% %but have sum of pii is 1
% pinv([(P-eye(n))'; ones(1,n)])*[zeros(n,1);1]
% 
% P^100 %another way



