Numits=30000;
Numsims=100;
window=200;

thetastar1 = 1.1830;
thetastar2=3.3918;

vH=50;
v=50;
weighty=2;
vp = v*weighty;

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);

ETp=zeros(1,Numits);
EPp=zeros(1,Numits);
ERp=zeros(1,Numits);
ETpav=zeros(1,Numits);
EPpav=zeros(1,Numits);
ERpav=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;
i1=0;
i2=0;
i3=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;
Qp=0;
Jp=0;
gammaoldp=gammamin;

for k=1:Numits
    %Get A = [Time Reward Power] 
    if k<=Numits/3
    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];
    elseif k<=Numits*2/3
    U1=rand(1);
    U2=rand(1);
    A=[1 0 0; 1+9*U1 min(2*(10*(U2+1)),rmax) 1+9*U1; 6 + 6*U1 10*U1*(U2+1) U1]; 
    else
    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];    
    end
    
    %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;
    
    
    
    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 = Prime
    %Row selection
    
    T=A(1,1);
    R=A(1,2);
    P=A(1,3);
    Y=weighty*(P-qav*T);
    wbest=-vp*R + Jp*T + Qp*Y;
    ibest=1;
    for i=2:3
        T=A(i,1);
        R=A(i,2);
        P=A(i,3);
        Y=weighty*(P-qav*T);
        w = -vp*R + Jp*T + Qp*Y;
        if (w<wbest)
            wbest=w;
            ibest=i;
        end
    end
    
    T=A(ibest,1);
    R=A(ibest,2);
    P=A(ibest,3);
    Y=weighty*(P-qav*T);
    
    
    
    ETp(k)=ETp(k)*(sims-1)/sims + T/sims;
    ERp(k)=ERp(k)*(sims-1)/sims + R/sims;
    EPp(k)=EPp(k)*(sims-1)/sims + P/sims;
    
    %Gamma
    gammap = gammaoldp + (vp*R-Jp*T-Qp*Y)/(gammaoldp*alpha*vp^2);
    if gammap > gammamax
        gammap = gammamax;
    elseif gammap < gammamin
        gammap = gammamin;
    end
    
    %VQ
    
    Qp = Qp + Y;
    if Qp>q*vp
        Qp=q*vp;
    elseif Qp<0
        Qp=0;
    end
    
    Jp = Jp + T - 1/gammap;
    if Jp<0
        Jp=0;
    end
    gammaoldp=gammap;
    
     
    %%%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);

ETpav(1)=ETp(1);
ERpav(1)=ERp(1);
EPpav(1)=EPp(1);

for k=2:Numits
   if (k<=window)
   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;
   
   ETpav(k)=ETpav(k-1)*(k-1)/k + ETp(k)/k;
   ERpav(k)=ERpav(k-1)*(k-1)/k + ERp(k)/k;
   EPpav(k)=EPpav(k-1)*(k-1)/k + EPp(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;
   else
       ETav(k)=ETav(k-1) - ET(k-window)/window + ET(k)/window;
       ERav(k)=ERav(k-1) - ER(k-window)/window + ER(k)/window;
       EPav(k)=EPav(k-1) - EP(k-window)/window + EP(k)/window;
       
       ETpav(k)=ETpav(k-1) - ETp(k-window)/window + ETp(k)/window;
       ERpav(k)=ERpav(k-1) - ERp(k-window)/window + ERp(k)/window;
       EPpav(k)=EPpav(k-1) - EPp(k-window)/window + EPp(k)/window;
       
       ETavgreedy(k)=ETavgreedy(k-1) - ETgreedy(k-window)/window + ETgreedy(k)/window;
       ERavgreedy(k)=ERavgreedy(k-1) - ERgreedy(k-window)/window + ERgreedy(k)/window;
       EPavgreedy(k)=EPavgreedy(k-1) - EPgreedy(k-window)/window + EPgreedy(k)/window;
       
       ETavH(k)=ETavH(k-1) - ETH(k-window)/window + ETH(k)/window;
       ERavH(k)=ERavH(k-1) - ERH(k-window)/window + ERH(k)/window;
       EPavH(k)=EPavH(k-1) - EPH(k-window)/window + EPH(k)/window;
   end
   
end

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

Rewardp = ERpav./ETpav;
Powerp = EPpav./ETpav;

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,Powerp, '-r', times, PowerH, 'b')
%ylim([0,1])
plot(times, Reward, '-g', times, RewardH, '-b', times, Rewardp, '-r', times, thetastar1*ones(size(times)), '--', times, thetastar2*ones(size(times)), '--')
ylim([0,5])


D2 = [Power;Reward;PowerH;RewardH;Powerp;Rewardp];
i1av = i1/(Numits/2)
i2av = i2/(Numits/2)
i3av=i3/(Numits/2)
i1av+i2av+i3av

i1avgreedy
i2avgreedy
i3avgreedy
i1avgreedy+i2avgreedy+i3avgreedy