% Figure 4.4 page 128 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 3.

clc
clear
v=0;
Vmin=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
];

T1=[0	0.875	0.0625	0.0625
0	1	0	0
1	0	0	0
];

T2=[0	0.75	0.125	0.125
0	1	0	0
1	0	0	0
];

T3=[0	0	0.5	0.5
0	1	0	0
1	0	0	0
];

T4=[0	0	0	1
0	1	0	0
1	0	0	0
];

Aph=0.9;
alph=0.2

%***************************************
%Step 0
V=zeros(4,1);
VV(1,:)=V(:,1)';
Policy=zeros(1,4);
ite=1;
S=1; % initial state step 0b.
%***************************************
%Step 1 and 2

if S==1
T=T1;
elseif S==2
    T=T2;
  elseif S==3
    T=T3;
    elseif S==4
    T=T4;
end

for a =1 :1:3
tempv(1,a)=C(S,a)+(Aph*(T(a,:)*V));
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 3

for state=1:1:4
if state==S
V(state,1)=(1-alph)*V(state,1)+alph*v
end
end

% alpha decay
ee=5000000*2;
    uz=ite^2/(ee+ite);
    alph=alph/(1+uz);

%step4

T11=[0	0.875	0.9375	1
0	1	1	1
1	1	1	1
];

T22=[0	0.75	0.825	1
0	1	1	1
1	1	1	1
];

T33=[0	0	0.5 1
0	1	1	1
1	1	1	1
];

T44=[0	0	0	1
0	1	1	1
1	1	1	1
];

if S==1
TT=T11;
elseif S==2
    TT=T22;
  elseif S==3
    TT=T33;
    elseif S==4
    TT=T44;
end

for a=1:1:3
    if a==Act
        temp= TT(Act,:);
    end
end

f=rand(1);
if f>temp(1,1) & f<=temp(1,2)
    Snext=2;
elseif f>temp(1,2) & f<=temp(1,3)
    Snext=3;
elseif f>temp(1,3) & f<=temp(1,4)
    Snext=4;
else
    Snext =1;
end
VV(2,:)=V(:,1)';
%***********************************************************
%Step 3
n=3;
while ite<500
ite=ite+1;
S=Snext;

%Step 1 and 2



if S==1
T=T1;
elseif S==2
    T=T2;
  elseif S==3
    T=T3;
    elseif S==4
    T=T4;
end

for a =1 :1:3
tempv(1,a)=C(S,a)+(Aph*(T(a,:)*V));
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 3

for state=1:1:4
if state==S
V(state,1)=(1-alph)*V(state,1)+alph*v
end
end

% alpha decay
ee=5000000*2;
    uz=ite^2/(ee+ite);
    alph=alph/(1+uz);
VV(n,:)=V(:,1)';
n=n+1;

%step4

T11=[0	0.875	0.9375	1
0	1	1	1
1	1	1	1
];

T22=[0	0.75	0.825	1
0	1	1	1
1	1	1	1
];

T33=[0	0	0.5 1
0	1	1	1
1	1	1	1
];

T44=[0	0	0	1
0	1	1	1
1	1	1	1
];

if S==1
TT=T11;
elseif S==2
    TT=T22;
  elseif S==3
    TT=T33;
    elseif S==4
    TT=T44;
end

for a=1:1:3
    if a==Act
        temp= TT(Act,:);
    end
end

f=rand(1);
if f>temp(1,1) & f<=temp(1,2)
    Snext=2;
elseif f>temp(1,2) & f<=temp(1,3)
    Snext=3;
elseif f>temp(1,3) & f<=temp(1,4)
    Snext=4;
else
    Snext =1;
end

end
%******************************************
Policy

plot(VV)