%This code is to study 4 types of alpha decay (too fast, fast, medium, and slow)
%and its relation to number of iterations
% The study is conducted by observing the number of times the Policy {1123}
% is generated and the average MSE

%1500
%ee=500000000*7; slow
%ee=50000000*7; medium
%ee=5000000*7; fast alpha died at 800
%1/ite is the too fast scheme

%500
%ee=50000000*2;slow
%ee=5000000*2; medium
%ee=500000*2; fast alpha died at 250
%1/ite is the too fast scheme

% Figure 4.7 page 141 of text 
% Instead of V bar (approximation), V from previous iterations are stored
% and used in place of V bar.

% The algorithm avoids looping over all states.
% Stopping criteria is the number of iterations, which is fixed to an arbitrarily large number.
% This is implemented in infinite horizon. Therefore there is no time t
% Alpha decay is used for step 2b.

% use of trasition probaility is avoided
% Post decision state variable is introduced
% v is calculated around post decision state variable

clc
clf
clear
% number of simulations
counter=1;
simnum=1;
Store_policy=0;
MSQEE=0;
Mean_MSQEE=0;
Std_MSQEE= 0;
Min_MSQEE=0;
Max_MSQEE=0;
alphs=0;
while simnum<=100
    
v=0;
Vmin=0;
msqe=0;
gradient =0;

%State 0,1,2,3 is 1,2,3,4 for Matlab
C=[	0	4000	6000
	1000	4000	6000
	3000	4000	6000
	10000	10000	6000
];

% post desion state PDS. This is a (pre-decision state X action) matrix
PDS= [1 2 1 
      2 2 1
      3 2 1
      4 2 1];

Aph=0.95;
alph=0.7;

%***************************************
%Step 0
V=zeros(4,1);
VV(1,:)=V(:,1)';
Policy=zeros(1,4);
ite=1;
S=1; % initial state
xSx= 1; % initial Post decision state
%***************************************
%Step 1 and 2

% Step 2a  

for a =1 :1:3
    Sx= PDS(S,a);
tempv(1,a)=C(S,a)+(Aph*V(Sx,1));
end
v=min(tempv(1,:));

for a =1 :1:3
if v-tempv(1,a)<0.001 & v-tempv(1,a)>-0.001
Act=a;
end
end

for state=1:1:4
if state==S
Policy(1,state)=Act;
end
end

%Step 2b



gradient(ite,1)=V(xSx,1)-v;
gradient(ite,2)=(V(xSx,1)-v)^2;
msqe(ite,1)= sum(gradient(:,2))/ite;

V(xSx,1)=(1-alph)*V(xSx,1)+alph*v;
xSx=PDS(S,Act); % Find the post decision state

%Step 2c
%Finding the predecision state (This is a random process)
f=rand(1);
if f<0.25
    Snext=1;
elseif f>=0.25 & f<0.5
    Snext =2;
elseif f>=0.5 & f<0.75
    Snext =3;
elseif f>=0.75
    Snext=4;

end
VV(2,:)=V(:,1)';
%***********************************************************
%Step 3
n=3;
while ite<1500
ite=ite+1;
S=Snext;

% alpha decay schemes
%1500
%ee=500000000*7; slow
%ee=50000000*7; medium
%ee=5000000*7; fast alpha died at 800
%1/ite is the too fast scheme

ee=5000000000*7;
    uz=ite^2/(ee+ite);
    alph=alph/(1+uz);   
      alphs(ite,1)=alph;
      
      
%Step 1 and 2

% Step 2a is implemented with the assumption that all outcomes are
% possible. Therefore summation in step 2b is over all possible states as
% given in figure 4.1 page 97.

%Step 2b
for a =1 :1:3
    Sx= PDS(S,a);
tempv(1,a)=C(S,a)+(Aph*V(Sx,1));
end
v=min(tempv(1,:));

for a =1 :1:3
if v-tempv(1,a)<0.001 & v-tempv(1,a)>-0.001
Act=a;
end
end

for state=1:1:4
if state==S
Policy(1,state)=Act;
end
end

%Step 2b


gradient(ite,1)=V(xSx,1)-v;
gradient(ite,2)=(V(xSx,1)-v)^2;
msqe(ite,1)= sum(gradient(:,2))/ite;

V(xSx,1)=(1-alph)*V(xSx,1)+alph*v;
xSx=PDS(S,Act); % Find the post decision state


VV(n,:)=V(:,1)';
n=n+1;
%Step 2c

%Finding the predecision state (This is a random process)
f=rand(1);
if f<0.25
    Snext=1;
elseif f>=0.25 & f<0.5
    Snext =2;
elseif f>=0.5 & f<0.75
    Snext =3;
elseif f>=0.75
    Snext=4;

end

end
Store_policy(simnum,1)=Policy(1,1);
Store_policy(simnum,2)=Policy(1,2);
Store_policy(simnum,3)=Policy(1,3);
Store_policy(simnum,4)=Policy(1,4);
if Policy(1,1) == 1;
    if  Policy(1,2) == 1;
         if Policy(1,3) == 2;
         if  Policy(1,4) == 3;
          counter=counter+1;
      end
  end
end
end
MSQEE(simnum,1)=msqe(ite,1);
simnum=simnum+1
end
%******************************************

% Plots for the last simulation

% subplot(3,1,1), plot(VV)
% subplot(3,1,2),plot (msqe)
% subplot(3,1,3),plot (MSQEE)
plot(VV)
figure
plot (msqe)
figure
plot(alphs)
Store_policy
counter
Mean_MSQEE=mean(MSQEE)
Std_MSQEE= std(MSQEE)
Min_MSQEE=min(MSQEE)
Max_MSQEE=max(MSQEE)