% Value interation for average cost dynamic program
N=200;  % the number of points in state space in the unit interval
M=400;  % the number of points in action space in the unit interval
Iteration=20; % the mumber of interation
z=linspace(0,1,N); 
dz=1/(N-1);
gama=linspace(0,1,M);
dgama=1/(M-1);
delta=linspace(0,1,M);
J=zeros(1,N);
delta_max=zeros(1,N);
gama_max=zeros(1,N);
for index_iteration=1:Iteration
    Jtemp=zeros(1,N);
    for index_z=2:N-1
        J_delta_gama=zeros(M,M);
        Index_gama=floor((1-z(index_z))/dgama)+1;     % gama is less than 1-z
        Index_delta=floor(z(index_z)/dgama)+1;        % delta is less than z
        for index_delta=1:Index_delta
            for index_gama=1:Index_gama
                d=delta(index_delta);
                r=gama(index_gama);
                q=1/2+(d-r)/2;
                Part1=-q*log2(q)-(1-q)*log2(1-q)+d+r-1;
                Index=round(2*d/(1+d-r)/dz);
                if Index<1;
                    Index=Index+1;
                end
                Part2=(1+d-r)/2*J(Index);
                Index=round((1-2*r/(1-d+r))/dz);
                if Index<1;
                    Index=Index+1;
                end
                Part3=(1-d+r)/2*J(Index);
                J_delta_gama(index_delta,index_gama)=Part1+Part2+Part3;
            end
        end
        [C,I]=max(J_delta_gama);
        [C2,I2]=max(C);
        Jtemp(index_z)=C2;
        Index_Max_delta=I(I2);
        Index_Max_gama=I2;
        delta_max(index_z)=(Index_Max_delta-1)*dgama;
        gama_max(index_z)=(Index_Max_gama-1)*dgama;
%         Jtemp(index_z)=max(max(J_delta_gama));
    end
    figure
    plot(delta_max)
    figure
    plot(gama_max)
    J=Jtemp;
    figure
    plot(J)
end