% 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=[	50	53.4	65.2	86.2
18.4	40.2	61.2	10000
5.2	36.2	10000	10000
1.2	10000	10000	10000
];

% post desion state PDS. This is a (pre-decision state on left and action
% on top) matrix, the contents of which are PDS.
PDS= [1 2 3 4
      2 3 4 4
      3 4 4 4
      4 4 4 4];

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 step 0c
xSx= 1; % initial Post decision state
%***************************************
%Step 1 and 2

% Step 2a  

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

% Identify the action that minimized v in the previous step
for a =1 :1:4
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


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


%Step 2c
%Finding the next predecision state with uniform demand prob
f=rand(1);
if f<0.25
    Snextt=0;
elseif f>=0.25 & f<0.5
    Snextt =1;
elseif f>=0.5 & f<0.75
    Snextt =2;
elseif f>=0.75
    Snextt=3;
end
Snext= max((xSx-(Snextt)), 1);
VV(2,:)=V(:,1)';
%***********************************************************
%Step 3
n=3;
while ite<1500
ite=ite+1;
S=Snext;
% alpha decay schemes

%6000
%5000000000*7
%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=50000000*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:4
    Sx= PDS(S,a);
tempv(1,a)=C(S,a)+(Aph*V(Sx,1));
end
v=min(tempv(1,:));

for a =1 :1:4
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 using uniform demand prob
f=rand(1);
if f<0.25
    Snextt=0;
elseif f>=0.25 & f<0.5
    Snextt =1;
elseif f>=0.5 & f<0.75
    Snextt =2;
elseif f>=0.75
    Snextt=3;

end
%Finding the predecision state (This is a random process)
Snext= max((xSx-(Snextt)), 1);
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) == 4;
    if  Policy(1,2) == 3;
         if Policy(1,3) == 1;
         if  Policy(1,4) == 1;
          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-1
Mean_MSQEE=mean(MSQEE)
Std_MSQEE= std(MSQEE)
Min_MSQEE=min(MSQEE)
Max_MSQEE=max(MSQEE)
