Numits=10000;
Numsims = 40;

v=10;


RobbinsdataT=zeros(1,Numits);
RobbinsdataR=zeros(1,Numits);
greedydataT=zeros(1,Numits);
greedydataR=zeros(1,Numits);
dataT=zeros(1,Numits);
dataR=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);

p0=.1;
p=0.3;
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;



D1RobbinsT=[];
D1greedyT=[];
D1T=[];
D1RobbinsR=[];
D1greedyR=[];
D1R=[];

J=0;
gammaold=gammamin; 
for k=1:Numits
    %Get project matrix A
    if (k<Numits/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
       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 
    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);
    
    
    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; 
    
D1RobbinsT=[D1RobbinsT TavRobbins];   
D1greedyT = [D1greedyT Tavgreedy];
D1T=[D1T Tav];
D1RobbinsR=[D1RobbinsR RavRobbins];   
D1greedyR = [D1greedyR Ravgreedy];
D1R=[D1R Rav];

ERobbinsdataT(k)=ERobbinsdataT(k)*(sims-1)/sims + D1RobbinsT(k)/sims;
EgreedydataT(k)=EgreedydataT(k)*(sims-1)/sims + D1greedyT(k)/sims;
EdataT(k)=EdataT(k)*(sims-1)/sims + D1T(k)/sims;
ERobbinsdataR(k)=ERobbinsdataR(k)*(sims-1)/sims + D1RobbinsR(k)/sims;
EgreedydataR(k)=EgreedydataR(k)*(sims-1)/sims + D1greedyR(k)/sims;
EdataR(k)=EdataR(k)*(sims-1)/sims + D1R(k)/sims;


end

end

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

eta
theta
performanceRobbins
performancegreedy
performanceadapt
TavRobbins
RavRobbins
J
gamma
Tav
Rav

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

ValdataRobbins=ERobbinsdata;
Valdatagreedy=Egreedydata;
