% Figure 4.7 page 141 of text 
%  V bar (approximation) is used (4 different schemes are evaluated)
%Section 8.2.3 on page 313 is used to update Theta cofficients

% 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
clear
clf
counter=0;
counter1=0;
simnum=1;
MSQEE=0;
Mean_MSQEE=0;
Std_MSQEE= 0;
Min_MSQEE=0;
Max_MSQEE=0;
Store_policy=0;
Sxx=0;
X=0;
while simnum<=100
v=0;
Vmin=0;
Theta=0;
Gama=0;
Vbar=0;
Sxx=0;
xSx=0;
t=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
	15000	15000	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];
 % basis function with delta = 2 (std dev of the states) and 
 %Data - [1
%   2
%   3
%   4
%   ];
   X= [     
        -0.7950    0.4812   -0.3165   -0.1904
   -0.5432   -0.3195    0.5480    0.5500
   -0.2566   -0.6558    0.0721   -0.7063
   -0.0838   -0.4860   -0.7710    0.4030

];


Aph=0.9;
alph=0.7;
NN=4; %number of states
%***************************************
%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;

% alpha decay
    ee=50000000*7;
    uz=ite^2/(ee+ite);
    alph=alph/(1+uz);
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

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

%Step 1 and 2

%Function approximation
% let phi(Sx)=Sx
% ensuring all states have been visited atleast once. This step is not
% needed when you sample states using exploration from a large number of states.

for p=1:1:NN
    if V(p,1)~=0
        if Ccount~=5
        Ccount=Ccount+1;
        end
    end
end
if Ccount<4
    Ccount=0;
end

if Ccount==4;
    %t=[1 2 3 4]'; %system states
  % X = [ones(size(t))  exp(-t) t.*exp(-t) ];
 % X = [ones(size(t))  exp(-t) ];
 % X = [ones(size(t))  t ];
 % X = [ones(size(t))  t t.*t];

Theta = X\V;
Ccount=5;
end

%Step 2a
for a =1 :1:3
    Sx= PDS(S,a);
 %   Sxx=[ 1 exp(-Sx) Sx.*exp(-Sx)];
%   Sxx=[ 1 exp(-Sx) ];
%   Sxx=[ 1 Sx];
%   Sxx=[ 1 Sx Sx^2];
gg=size(X);
g=gg(1,2);
for h=1:1:g
Sxx(1,h)=X(Sx,h);
end
    if Ccount==5;
    Vbar = ((Theta')/2) * Sxx' ;
        V(Sx,1)=Vbar;
    end
    
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
 %    Sxx=[ 1 exp(-xSx) Sx.*exp(-xSx)];
 %   Sxx=[ 1 exp(-xSx) ];
 %   Sxx=[ 1 xSx];
 %   Sxx=[ 1 xSx xSx^2];
for h=1:1:g
Sxx(1,h)=X(xSx,h);
end
if Ccount==5;
    Vbar = ((Theta')/2) * Sxx' ;
    V(xSx,1)=Vbar;
end
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;

%theta update
% if Ccount==5;
% Theta= Theta- alph*(V(xSx,1)-v)*Sxx';
% end
Theta = X\V;

% alpha decay
    ee=50000000*7;
    uz=ite^2/(ee+ite);
    alph=alph/(1+uz);
VV(n,:)=V(:,1)';
n=n+1;
xSx=PDS(S,Act); % Find the next 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
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
if Policy(1,1) == 1;
    if  Policy(1,2) == 1;
         if Policy(1,3) == 1;
         if  Policy(1,4) == 3;
          counter1=counter1+1;
      end
  end
end
end
MSQEE(simnum,1)=msqe(ite,1);
simnum=simnum+1
end
%******************************************
%subplot(3,1,1), plot(VV)
subplot(2,1,1),plot (msqe)
subplot(2,1,2),plot (MSQEE)
Store_policy
counter
counter1
Mean_MSQEE=mean(MSQEE)
Std_MSQEE= std(MSQEE)
Min_MSQEE=min(MSQEE)
Max_MSQEE=max(MSQEE)


