function [h_optt] = Theoritical_opt(lambda,B,c_f,c_delta)
%% Theoritical optimization 


%% Iterative identification of active constraints
N=length(lambda)
flag=1;
h_optt=zeros(N,1); %optimal hit prob
differenc_set=1:N; % the set of variables to optimize 
fix_i=[]; % set of variables that reached optimality and are excluded from the optimization
n=N; 
b=B;

%%
while flag
    if(n==0) 
        if(b>0) % if there is left over cache size and mu is not zero (the loop would break), redistribute it among the zero hit probability
            differenc_set=find(h_optt==0)';
            fix_i=setdiff(1:N,differenc_set)';
            n=length(differenc_set);
            continue;     
        else           
             h_optt(differenc_set)=0;
             break;
        end
    end
    % Optimal Lagrangian mult. and hit prob. calculated theoritically for the set of variables in  differenc_set
    mu=max(0,(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set)));  %optimal lagrangian mult.
    h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta %optimal hit prob
    
    % mu has to be >=0
    if(mu<0)
        b=(n*c_f/c_delta); % this sets mu to zero in the next iteration
        continue;
    end
    
    % check the violation of the hit_prob const  
    larger_i=find(h_optt>1); % h>1
    smaller_i=find(h_optt<0); % h<0

%     smaller=h(find(differenc_set<0))-0;
    % no violation means optimal solution is reached for all objects
    if(length(smaller_i)==0 && length(larger_i)==0)
%         flag=0;
  
        break;
    end
    
    % find the furthest object from the 0 boundary
    min_viol=0;
    min_viol_i=-1;
    if(length(smaller_i)>0)
        [min_viol, min_viol_i]=min(h_optt);
    end
     % find the furthest object from the 1 boundary
    max_viol=0; 
    max_viol_i=-1;
    if(length(larger_i)>0)
       larger=h_optt-1;
       [max_viol ,max_viol_i]=max(h_optt-1);

    end
    
     % compare both furthest objects from both boundaries
    viol_i=min_viol_i;
    min_viol_flag=1; % True if the furthest one is from the left
    if(max_viol>abs(min_viol))
        viol_i= max_viol_i;
        min_viol_flag=0;

    end
     % set the furthest object to the nearest boundary
     if(min_viol_flag)
        h_optt(viol_i)=0;
    else
        h_optt(viol_i)=min(1,b);
     end
    
    %calculate the new parameters after removing the furthest object from
    %the decision variables
    B_new=b-(h_optt(viol_i));
    b=B_new;
    fix_i=[fix_i' viol_i']';
    differenc_set=setdiff(1:N,fix_i)   ; 
    n=N-length(fix_i);
    
    
    %     % Identify the most violating object from the right side  h>1
%     if(length(larger_i)>0)
%        larger=h_optt-1;
%        [max_viol ,max_viol_i]=max(h_optt-1); % maximum violating object
%        h_optt(max_viol_i)=min(1,b); %project to the feasible range
%        b=max(b-1,0); % update the cache size 
%        fix_i=[fix_i' max_viol_i']'; %exclude i from the set of decision variables
%        differenc_set=setdiff(1:N,fix_i); % obtain the set of decision variables 
%        n=N-length(fix_i); % update the number of decision variables 
%        continue;
%     end
%     
%     if(length(smaller_i)>0)
%         [min_viol, min_viol_i]=min(h_optt);
%         h_optt(min_viol_i)=0; 
%         fix_i=[fix_i' min_viol_i']';
%         differenc_set=setdiff(1:N,fix_i)   ; 
%         n=N-length(fix_i);
%     end
%     
% end

end

%% Identfying the active constraints collectively
% flag=1;
% h_optt=zeros(N,1);
% differenc_set=1:N;
% fix_i=[];
% n=N;
% b=B;
% while flag
%     mu=(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set));
%     h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta
% 
%     larger_i=find(h_optt>1);    
% %     larger=h_optt(larger_i)-1;
%     smaller_i=find(h_optt<0);
% %     smaller=h(find(differenc_set<0))-0;
%     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
% 
%     if(length(larger_i)+length(smaller_i)==0)
%         flag=0;
%         break;
%     end
%     if(length(smaller_i)>0)
%         h_optt(smaller_i)=0;
%         fix_i=[fix_i' smaller_i' ]';
%         differenc_set=setdiff(1:N,fix_i)    
%         n=N-length(fix_i);
%         continue
%     end
% %     h_optt(smaller_i)=0;
%     if(length(larger_i)>b)
%         [~,index]=maxk(h_optt,b)
%         h_optt(index)=1;
%         B_new=b-sum(h_optt(index));
%         fix_i=[fix_i' smaller_i' index']';
%     else
%          h_optt(larger_i)=1;
%          B_new=b-sum(h_optt(larger_i));
%          fix_i=[fix_i' smaller_i' larger_i']';
%     end
% %     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
% 
%     b=B_new;
%    
%     differenc_set=setdiff(1:N,fix_i)    
%     n=N-length(fix_i);
% end
% % h_optt=zeros(N,1);
% % h_optt(end-B+1:end)=1;
% optimal=[h_opt h_optt]
% obj_1=sum(lambda .* ((1-h_opt)*c_f+h_opt.^2*c_delta/2));
% obj_2=sum(lambda .* ((1-h_optt)*c_f+h_optt.^2*c_delta/2));
% objective=[obj_1 obj_2]
% const_1=sum(h_opt)-B;
% const_2=sum(h_optt)-B;
% constraint=[const_1 const_2]
% 


end

