Numits=20000;
Numsims = 40;

star1=33.7469;
star2 = 54.6846;

v=10;
window=200;

RobbinsdataT=zeros(1,Numits);
RobbinsdataR=zeros(1,Numits);
greedydataT=zeros(1,Numits);
greedydataR=zeros(1,Numits);
dataT=zeros(1,Numits);
dataR=zeros(1,Numits);

ERobbinsdata=zeros(1,Numits);
Egreedydata=zeros(1,Numits);
Edata=zeros(1,Numits);
ERobbinsdataT=zeros(1,Numits);
EgreedydataT=zeros(1,Numits);
EdataT=zeros(1,Numits);
ERobbinsdataR=zeros(1,Numits);
EgreedydataR=zeros(1,Numits);
EdataR=zeros(1,Numits);


rmax=50*10;
tmin=1;
tmax = 10;
gammamin=1/tmax;
gammamax=1/tmin;
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

thetamin = 0;
thetamax = 500;

%%%%
for sims=1:Numsims

Tvec=zeros(1,window);
Rvec=zeros(1,window);
Tgreedyvec=zeros(1,window);
Rgreedyvec=zeros(1,window);
TRobbinsvec=zeros(1,window);
RRobbinsvec=zeros(1,window);
    
theta=0;
TavRobbins=0;
RavRobbins=0;
Tavgreedy=0;
Ravgreedy=0;
Tav = 0 ;
Rav = 0;



D1Robbins=[];
D1greedy=[];
D1=[];

J=0;
gammaold=gammamin; 
for k=1:Numits
    %Get project matrix A
    if (k<Numits/2)
        p0=.1;
        p=0.3;
        p1=1-p0-p;
        p2 = p/2;
        p3=p/2;
    A = [1 0];
    X=rand();
    if X<p0
        M=1;
    elseif X<p0+p1
        M=2;
    elseif X<p0+p1+p2
        M=3;
    else
        M=4;
    end
    for j=1:(M-1)
       X = rand()*50; 
       T = 1.0+rand()*9;
       A = [A; [T X*T]];
    end
    else
        p0=0;
        p=0.8;
        p1=1-p0-p;
        p2 = p/2;
        p3=p/2;
       A = [1 0];
    X=rand();
    if X<p0
        M=1;
    elseif X<p0+p1
        M=2;
    elseif X<p0+p1+p2
        M=3;
    else
        M=4;
    end
    for j=1:(M-1)
       X = 10+rand()*20; 
       Y = rand()*200;
       T = 1.0+rand()*9;
       A = [A; [T (X*T+Y)]];
    end 
    end
    
    %Robbins
    bestj = 1;
    bestweight = A(1,2)-theta*A(1,1);
    for j=2:M
        weight = A(j,2)-theta*A(j,1);
        if weight > bestweight
            bestweight=weight;
            bestj=j;
        end
    end
    
    TRobbins = A(bestj,1);
    RRobbins = A(bestj,2);
    %%%window
    if (k<=window)
       TRobbinsvec(k)=TRobbins;
       RRobbinsvec(k)=RRobbins;
       if (k==1)
           D2RobbinsT(k)=TRobbins;
           D2RobbinsR(k)=RRobbins;
       else
          D2RobbinsT(k)=D2RobbinsT(k-1)*(k-1)/k+TRobbins/k; 
          D2RobbinsR(k)=D2RobbinsR(k-1)*(k-1)/k+RRobbins/k;
       end
    else
        D2RobbinsT(k)=D2RobbinsT(k-1) - TRobbinsvec(1)/window+TRobbins/window;
        D2RobbinsR(k)=D2RobbinsR(k-1) - RRobbinsvec(1)/window+RRobbins/window; 
        TRobbinsvec = [TRobbinsvec(2:window) TRobbins];
        RRobbinsvec = [RRobbinsvec(2:window) RRobbins];
    end
    %%%
    
    
    eta = 1/(k+1);
    theta = theta + eta*(RRobbins-theta*TRobbins); 
    if theta > thetamax
        theta = thetamax;
    elseif theta < thetamin
        theta = thetamin;
    end
    
    TavRobbins = TavRobbins*(k-1)/k + TRobbins*1/k;
    RavRobbins = RavRobbins*(k-1)/k + RRobbins*1/k; 
    
    %Greedy
    bestj = 1;
    bestweight = A(1,2)/A(1,1);
    for j=2:M
        weight = A(j,2)/A(j,1);
        if weight > bestweight
            bestweight=weight;
            bestj=j;
        end
    end
    
    Tgreedy = A(bestj,1);
    Rgreedy = A(bestj,2);
    
    
    Tavgreedy = Tavgreedy*(k-1)/k + Tgreedy*1/k;
    Ravgreedy = Ravgreedy*(k-1)/k + Rgreedy*1/k; 
    
    %%%window
    if (k<=window)
       Tgreedyvec(k)=Tgreedy;
       Rgreedyvec(k)=Rgreedy;
       if (k==1)
           D2greedyT(k)=Tgreedy;
           D2greedyR(k)=Rgreedy;
       else
          D2greedyT(k)=D2greedyT(k-1)*(k-1)/k+Tgreedy/k; 
          D2greedyR(k)=D2greedyR(k-1)*(k-1)/k+Rgreedy/k;
       end
    else
        D2greedyT(k)=D2greedyT(k-1) - Tgreedyvec(1)/window+Tgreedy/window;
        D2greedyR(k)=D2greedyR(k-1) - Rgreedyvec(1)/window+Rgreedy/window; 
        Tgreedyvec = [Tgreedyvec(2:window) Tgreedy];
        Rgreedyvec = [Rgreedyvec(2:window) Rgreedy];
    end
    %%%
    
    %adaptive
    bestj = 1;
    bestweight = v*A(1,2)-J*A(1,1);
    for j=2:M
        weight = v*A(j,2)-J*A(j,1);
        if weight > bestweight
            bestweight=weight;
            bestj=j;
        end
    end
    T=A(bestj,1);
    R=A(bestj,2);
    Tav = Tav*(k-1)/k + T*1/k;
    Rav = Rav*(k-1)/k + R*1/k;
    
    gamma = gammaold + (v*R - J*T)/(gammaold*alpha*v^2);
    if gamma>gammamax
        gamma = gammamax;
    elseif gamma < gammamin
        gamma = gammamin;
    end
    
    J = J + T - 1/gamma;
    if J<0
        J=0;
    end
    gammaold = gamma; 
    
     %%%window
    if (k<=window)
       Tvec(k)=T;
       Rvec(k)=R;
       if (k==1)
           D2T(k)=T;
           D2R(k)=R;
       else
          D2T(k)=D2T(k-1)*(k-1)/k+T/k; 
          D2R(k)=D2R(k-1)*(k-1)/k+R/k;
       end
    else
        D2T(k)=D2T(k-1) - Tvec(1)/window+T/window;
        D2R(k)=D2R(k-1) - Rvec(1)/window+R/window; 
        Tvec = [Tvec(2:window) T];
        Rvec = [Rvec(2:window) R];
    end
    %%%
    
    

ERobbinsdataT(k)=ERobbinsdataT(k)*(sims-1)/sims + D2RobbinsT(k)/sims;
ERobbinsdataR(k)=ERobbinsdataR(k)*(sims-1)/sims + D2RobbinsR(k)/sims;

EgreedydataT(k)=EgreedydataT(k)*(sims-1)/sims + D2greedyT(k)/sims;
EgreedydataR(k)=EgreedydataR(k)*(sims-1)/sims + D2greedyR(k)/sims;

EdataT(k)=EdataT(k)*(sims-1)/sims + D2T(k)/sims;
EdataR(k)=EdataR(k)*(sims-1)/sims + D2R(k)/sims;



end

end

Egreedydata = EgreedydataR./EgreedydataT;
ERobbinsdata = ERobbinsdataR./ERobbinsdataT;
Edata = EdataR./EdataT;

performanceRobbins=RavRobbins/TavRobbins;
performancegreedy=Ravgreedy/Tavgreedy;
performanceadapt=Rav/Tav;

eta
theta
performanceRobbins
performancegreedy
performanceadapt
TavRobbins
RavRobbins
J
gamma
Tav
Rav


%timesnew=[49000:60000];
%ERobbinsdatanew=ERobbinsdata(49000:60000);
%Egreedydatanew=Egreedydata(49000:60000);
%Edatanew=Edata(49000:60000);
%tstar = ERobbinsdata(Numits);
%times = [1:Numits];
%plot(timesnew,ERobbinsdatanew, 'b', timesnew, Egreedydatanew, timesnew, Edatanew, '-g',timesnew, tstar*ones(size(Edatanew)),'--')


tstar = ERobbinsdata(Numits);
times = [1:Numits];
times1=[1:Numits/2];
times2=[Numits/2:Numits];
plot(times,ERobbinsdata, 'b', times, Egreedydata, '-r', times, Edata, '-g',times, star1*ones(size(times)),'--', times, star2*ones(size(times)), '--')