Numits=5000;
Numsims=40;
window=200;

thetastar = 1.1830;

vH=50;
v=200;

tmin=1;
tmax=12;
rmax=20;
pmax = 1;
pmin=0;
qav=1/3;
ymin = (-1)*(pmin - qav*tmax);
ymax=pmax-qav*tmin; 

gammamin=1/tmax;
gammamax=1/tmin;

q=100;
c1=rmax + (tmax-tmin)*(1+rmax)/tmin;
c2=((tmax-tmin)/tmin)*(tmax/tmin + tmin/tmax-2);
if c2>1/2
    alpha = c1/c2;
else
    alpha = 2*c1;
end


ET=zeros(1,Numits);
EP=zeros(1,Numits);
ER=zeros(1,Numits);
ETav=zeros(1,Numits);
EPav=zeros(1,Numits);
ERav=zeros(1,Numits);

ETgreedy=zeros(1,Numits);
EPgreedy=zeros(1,Numits);
ERgreedy=zeros(1,Numits);
ETavgreedy=zeros(1,Numits);
EPavgreedy=zeros(1,Numits);
ERavgreedy=zeros(1,Numits);


ETH=zeros(1,Numits);
EPH=zeros(1,Numits);
ERH=zeros(1,Numits);
ETavH=zeros(1,Numits);
EPavH=zeros(1,Numits);
ERavH=zeros(1,Numits);


for sims=1:Numsims
gammaold=gammamin;
Q=0;
J=0;
i1av=0;
i2av=0;
i3av=0;
i1avgreedy=0;
i2avgreedy=0;
i3avgreedy=0;
i1Hav=0;
i2Hav=0;
i3Hav=0;

QH=0;
theta=0;
TavH=0;
RavH=0;
for k=1:Numits
    %Get A = [Time Reward Power] 
    U1=rand(1);
    U2=rand(1);
    A=[1 0 0; 1+9*U1 10*U1*(U2+1) 1+9*U1; 6 + 6*U1 10*U1*(U2+1) U1];
    
    %Row selection
    
    T=A(1,1);
    R=A(1,2);
    P=A(1,3);
    Y=P-qav*T;
    wbest=-v*R + J*T + Q*Y;
    ibest=1;
    for i=2:3
        T=A(i,1);
        R=A(i,2);
        P=A(i,3);
        Y=P-qav*T;
        w = -v*R + J*T + Q*Y;
        if (w<wbest)
            wbest=w;
            ibest=i;
        end
    end
    
    T=A(ibest,1);
    R=A(ibest,2);
    P=A(ibest,3);
    Y=P-qav*T;
    
    if ibest==1
       i1av = i1av*(k-1)/k + 1/k;
    else
       i1av = i1av*(k-1)/k;
    end
    if ibest==2
       i2av = i2av*(k-1)/k + 1/k;
    else
       i2av = i2av*(k-1)/k;
    end
    if ibest==3
       i3av = i3av*(k-1)/k + 1/k;
    else
       i3av = i3av*(k-1)/k;
    end
    
    ET(k)=ET(k)*(sims-1)/sims + T/sims;
    ER(k)=ER(k)*(sims-1)/sims + R/sims;
    EP(k)=EP(k)*(sims-1)/sims + P/sims;
    
    %Gamma
    gamma = gammaold + (v*R-J*T-Q*Y)/(gammaold*alpha*v^2);
    if gamma > gammamax
        gamma = gammamax;
    elseif gamma < gammamin
        gamma = gammamin;
    end
    
    %VQ
    
    Q = Q + Y;
    if Q>q*v
        Q=q*v;
    elseif Q<0
        Q=0;
    end
    
    J = J + T - 1/gamma;
    if J<0
        J=0;
    end
    gammaold=gamma;
      
    %%% GREEDY
     %Row selection
    
    T=A(1,1);
    R=A(1,2);
    P=A(1,3);
    wbest=0;
    ibest=1;
    for i=2:3
        T=A(i,1);
        R=A(i,2);
        P=A(i,3);
        
        if (P/T<=qav)&(R/T>wbest)
            wbest=R/T;
            ibest=i;
        end
    end
    
    T=A(ibest,1);
    R=A(ibest,2);
    P=A(ibest,3);
    
    if ibest==1
       i1avgreedy = i1avgreedy*(k-1)/k + 1/k;
    else
       i1avgreedy = i1avgreedy*(k-1)/k;
    end
    if ibest==2
       i2avgreedy = i2avgreedy*(k-1)/k + 1/k;
    else
       i2avgreedy = i2avgreedy*(k-1)/k;
    end
    if ibest==3
       i3avgreedy = i3avgreedy*(k-1)/k + 1/k;
    else
       i3avgreedy = i3avgreedy*(k-1)/k;
    end
    
    ETgreedy(k)=ETgreedy(k)*(sims-1)/sims + T/sims;
    ERgreedy(k)=ERgreedy(k)*(sims-1)/sims + R/sims;
    EPgreedy(k)=EPgreedy(k)*(sims-1)/sims + P/sims;
    
    
    %%%H1
     %Row selection
    
    T=A(1,1);
    R=A(1,2);
    P=A(1,3);
    Y=P-qav*T;
    wbest=-vH*(R-theta*T)+ QH*Y;
    ibest=1;
    for i=2:3
        T=A(i,1);
        R=A(i,2);
        P=A(i,3);
        Y=P-qav*T;
        w = -vH*(R-theta*T)+ QH*Y;
        if (w<wbest)
            wbest=w;
            ibest=i;
        end
    end
    
    T=A(ibest,1);
    R=A(ibest,2);
    P=A(ibest,3);
    Y=P-qav*T;
    
    if ibest==1
       i1Hav = i1Hav*(k-1)/k + 1/k;
    else
       i1Hav = i1Hav*(k-1)/k;
    end
    if ibest==2
       i2Hav = i2Hav*(k-1)/k + 1/k;
    else
       i2Hav = i2Hav*(k-1)/k;
    end
    if ibest==3
       i3Hav = i3Hav*(k-1)/k + 1/k;
    else
       i3Hav = i3Hav*(k-1)/k;
    end
    
    ETH(k)=ETH(k)*(sims-1)/sims + T/sims;
    ERH(k)=ERH(k)*(sims-1)/sims + R/sims;
    EPH(k)=EPH(k)*(sims-1)/sims + P/sims;
    
        
    QH = QH + Y;
    if QH<0
        QH=0;
    end
    RavH = RavH*(k-1)/k + R/k;
    TavH = TavH*(k-1)/k + T/k;
    theta = RavH/TavH;
    
end
end

ETav(1)=ET(1);
ERav(1)=ER(1);
EPav(1)=EP(1);

ETavgreedy(1)=ETgreedy(1);
ERavgreedy(1)=ERgreedy(1);
EPavgreedy(1)=EPgreedy(1);

ETavH(1)=ETH(1);
ERavH(1)=ERH(1);
EPavH(1)=EPH(1);

for k=2:Numits
   ETav(k)=ETav(k-1)*(k-1)/k + ET(k)/k;
   ERav(k)=ERav(k-1)*(k-1)/k + ER(k)/k;
   EPav(k)=EPav(k-1)*(k-1)/k + EP(k)/k;
   
   ETavgreedy(k)=ETavgreedy(k-1)*(k-1)/k + ETgreedy(k)/k;
   ERavgreedy(k)=ERavgreedy(k-1)*(k-1)/k + ERgreedy(k)/k;
   EPavgreedy(k)=EPavgreedy(k-1)*(k-1)/k + EPgreedy(k)/k;
   
   ETavH(k)=ETavH(k-1)*(k-1)/k + ETH(k)/k;
   ERavH(k)=ERavH(k-1)*(k-1)/k + ERH(k)/k;
   EPavH(k)=EPavH(k-1)*(k-1)/k + EPH(k)/k;
   
end

Reward = ERav./ETav;
Power = EPav./ETav;

Rewardgreedy = ERavgreedy./ETavgreedy;
Powergreedy = EPavgreedy./ETavgreedy;

RewardH = ERavH./ETavH;
PowerH = EPavH./ETavH;

times=[1:Numits];
%plot(times, Power, '-g', times, qav*ones(size(times)), '--', times,Powergreedy, '-r', times, PowerH, 'b')
%ylim([0,1])
plot(times, Reward, '-g',times, Rewardgreedy, '-r', times, RewardH, '-b', times, thetastar*ones(size(times)), '--')
ylim([0,1.5])


D2 = [Power;Reward;PowerH;RewardH;Powergreedy;Rewardgreedy];
i1av
i2av
i3av
i1av+i2av+i3av

i1avgreedy
i2avgreedy
i3avgreedy
i1avgreedy+i2avgreedy+i3avgreedy