Numits=10000;
p0=.1;
p=.3;
Numsims = 40;

v=40;


window = 100;
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);


p1=1-p0-p;
p2 = p/2;
p3=p/2;
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

theta=0;
TavRobbins=0;
RavRobbins=0;
Tavgreedy=0;
Ravgreedy=0;
Tav = 0 ;
Rav = 0;

TRobbinsvec=zeros(1,window);
RRobbinsvec=zeros(1,window);
Tgreedyvec=zeros(1,window);
Rgreedyvec=zeros(1,window);
Tvec = zeros(1,window);
Rvec = zeros(1,window);

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

J=0;
gammaold=gammamin; 
for k=1:Numits
    %Get project matrix A
    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
    
    %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
           RobbinsdataT(k)=TRobbins;
           RobbinsdataR(k)=RRobbins;
        else
           RobbinsdataT(k)=RobbinsdataT(k-1)*(k-1)/k + TRobbins/k;
           RobbinsdataR(k)=RobbinsdataR(k-1)*(k-1)/k + RRobbins/k;
        end
    else
        RobbinsdataT(k) = RobbinsdataT(k-1) - TRobbinsvec(1)/window + TRobbins/window;
        RobbinsdataR(k) = RobbinsdataR(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; 
    
    %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; 
    
D1Robbins=[D1Robbins RavRobbins/TavRobbins];   
D1greedy = [D1greedy Ravgreedy/Tavgreedy];
D1=[D1 Rav/Tav];

ERobbinsdata(k)=ERobbinsdata(k)*(sims-1)/sims + D1Robbins(k)/sims;
Egreedydata(k)=Egreedydata(k)*(sims-1)/sims + D1greedy(k)/sims;
Edata(k)=Edata(k)*(sims-1)/sims + D1(k)/sims;


end

end

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

eta
theta
performanceRobbins
performancegreedy
performanceadapt
TavRobbins
RavRobbins
J
gamma
Tav
Rav

tstar = ERobbinsdata(Numits);
times = [1:Numits];
plot(times,ERobbinsdata, times, Egreedydata, times, Edata, times, tstar*ones(1,Numits),'--')