% prepare an upper bound of number of iterations for the algorithm
numiterssg_2_ub = 10000;
% initialize the counting on the actual number of iterations in the algorithm
step_count_ssg_2 = 0;

% prepare an 1*(numiterssg_2_ub+1) zero array to store the values of objective in each iteration
SSG_2_obj_candidate = zeros(1, numiterssg_2_ub+1);
% prepare an 1*(numiterssg_2_ub+1) zero array to store the values of s in each iteration
SSG_2_s_candidate = zeros(1, numiterssg_2_ub+1);
% prepare an 1*(numiterssg_2_ub+1) zero array to store the values of constraint in each iteration
SSG_2_cons_candidate = zeros(1, numiterssg_2_ub+1);
% prepare a d*(numiterssg_2_ub+1) zero matrix to store the iterates in each iteration by column
SSG_2_w_candidate = zeros(d, numiterssg_2_ub+1);

% prepare a d*(numiterssg_2_ub) zero matrix to store the subgradients of objective in each iteration by column
SSG_2_obj_subgrad_candidate = zeros(d, numiterssg_2_ub);
% prepare a d*(numiterssg_2_ub) zero matrix to store the subgradients of constraint in each iteration by column
SSG_2_cons_subgrad_candidate = zeros(d, numiterssg_2_ub);

% prepare an 1*(numiterssg_2_ub) zero array to store the step-sizes in each iteration
SSG_2_eta_candidate = zeros(1, numiterssg_2_ub);

% initialize I in SSG method
SSG_2_index_set = [];
% initialize CPU time counting in SSG method
SSG_2_time_set = [0];

% w_init is inherited from the hingeloss minimization, which is d*1 double
tic; % for recording CPU runtime
w = w_init;
SSG_2_w_candidate(:,1) = w_init; 
SSG_2_time_set(length(SSG_2_time_set)) = SSG_2_time_set(length(SSG_2_time_set)) + toc; % record CPU runtime

% change S here
S = -0.5:0.005:1.5; 
% prepare scaled S based on all scores by w = w^(0)
% this scaled S is used for all w^(t)'s
Sc_w = min(Ac*w) + ( max(Ac*w) - min(Ac*w) ) * S; 

% compute mean_male = mean( sigmoid ( Ac_male * w - Sc_w ) )
sigmoid_mean_male = mean( 1 ./ ( 1 + exp( Sc_w - Ac(idmalec,:)*w ) ) );
% compute mean_female = mean( sigmoid ( A_female * w - S_w ) )
sigmoid_mean_female = mean( 1 ./ ( 1 + exp( Sc_w - Ac(idfemalec,:)*w ) ) );
% get the array of objective candidates abs( mean_male - mean_female ) 
% then obtain the largest one with index
[obj, index] = max( abs( sigmoid_mean_male - sigmoid_mean_female ) ); 

% get the value of s
SSG_2_s_candidate(1,1) = S(1,index); 
% get the value of objective F(w^(0))
SSG_2_obj_candidate(1,1) = obj;

tic; % for recording CPU runtime
% compute the value of constraint G(w^(0))
SSG_2_cons_candidate(1,1) = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n - C;
SSG_2_time_set(length(SSG_2_time_set)) = SSG_2_time_set(length(SSG_2_time_set)) + toc; % record CPU runtime

fprintf('Iter=%i, s=%f, obj=%f, cons=%f\n', 0, SSG_2_s_candidate(1,1), SSG_2_obj_candidate(1,1)+0, SSG_2_cons_candidate(1,1)+0);
% add +0 to avoid the error: using fprintf. Function is not defined for sparse inputs

% define constants in the method
% Lambda = ( M*D_X + rho_hat*D_X^2 ) / C;
% epsilon = 0.001;
% E = 0.05;
% F = 3 * ( obj - 0 + rho_hat * D_X^2 / 2 ) / ( rho_hat * E ) + 1.5 * E * M^2;
E1 = 2e-4; % select from {5e-5, 1e-4, 2e-4, 5e-4}
E2 = 0.05; % select from {0.02, 0.05, 0.1, 0.2}

for t = 1:numiterssg_2_ub

    tic; % for recording CPU runtime
   
    ssg_2_epsilon_t = E1 / sqrt(t); 
    
    if SSG_2_cons_candidate(1,t)+0 < ssg_2_epsilon_t % perform subgradient on objective function

        % compute objevtive subgradient by calling the function "objective_subgradient"
        obj_subgrad = objective_subgradient(Ac, idmalec, idfemalec, nummalec, numfemalec, w, Sc_w);
        SSG_2_obj_subgrad_candidate(:,t) = obj_subgrad;
 
        ssg_2_eta_t =  E2 / sqrt(t); % compute eta_t
        SSG_2_eta_candidate(1,t) = ssg_2_eta_t;
        
        w = w - ssg_2_eta_t * obj_subgrad; % update w^(t)
        w = min( 0.5*D_X / norm(w,2), 1 ) * w; % projection to the ball with radius as 0.5*D_X
        SSG_2_w_candidate(:,t+1) = w;
        
        if t >= 1 + numiterssg_2_ub / 2
            SSG_2_index_set = [SSG_2_index_set, t-1]; % add t to I
        end

        SSG_2_time_set = [SSG_2_time_set, SSG_2_time_set(length(SSG_2_time_set))+toc]; % record CPU runtime

        step_count_ssg_2 = step_count_ssg_2 + 1;
        
    else % perform subgradient on constraint function

        % compute the constraint subgradient as the same as that in hingeloss_minimization
        cons_subgrad = - ( (A').*b' ) * ( 1 - b.*(A*w) > 0 ) / n;
        SSG_2_cons_subgrad_candidate(:,t) = cons_subgrad;

        ssg_2_eta_t = E2 / sqrt(t); % compute eta_t
        SSG_2_eta_candidate(1,t) = ssg_2_eta_t;
        
        w = w - ssg_2_eta_t * cons_subgrad; % update w^(t)
        w = min( 0.5*D_X / norm(w,2), 1 ) * w; % projection to the ball with radius as 0.5*D_X
        SSG_2_w_candidate(:,t+1) = w;
        
        SSG_2_time_set = [SSG_2_time_set, SSG_2_time_set(length(SSG_2_time_set))+toc]; % record CPU runtime

        step_count_ssg_2 = step_count_ssg_2 + 1;
        
    end

    % compute mean_male = mean( sigmoid ( Ac_male * w - Sc_w ) )
    sigmoid_mean_male = mean( 1 ./ ( 1 + exp( Sc_w - Ac(idmalec,:)*w ) ) );
    % compute mean_female = mean( sigmoid ( Ac_female * w - S_w ) )
    sigmoid_mean_female = mean( 1 ./ ( 1 + exp( Sc_w - Ac(idfemalec,:)*w ) ) );
    % get the array of objective candidates abs( mean_male - mean_female ) 
    % then obtain the largest one with index
    [obj, index] = max( abs( sigmoid_mean_male - sigmoid_mean_female ) ); 

    % get the value of s
    SSG_2_s_candidate(1,t+1) = S(1,index); 
    % get the value of objective F(w^(t))
    SSG_2_obj_candidate(1,t+1) = obj; 

    tic; % for recording CPU runtime
    % compute the value of constraint G(w^(t))
    SSG_2_cons_candidate(1,t+1) = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n - C;
    SSG_2_time_set(length(SSG_2_time_set)) = SSG_2_time_set(length(SSG_2_time_set)) + toc; % record CPU runtime

    fprintf('Iter=%i, s=%f, obj=%f, cons=%f\n', t, SSG_2_s_candidate(1,t+1), SSG_2_obj_candidate(1,t+1)+0, SSG_2_cons_candidate(1,t+1)+0);
    
end

% Option I: return w^(t) with the largest t in I
index_ssg_2_option_I = SSG_2_index_set(1,length(SSG_2_index_set));
s_ssg_2_option_I = SSG_2_s_candidate(1,index_ssg_2_option_I+1);
obj_ssg_2_option_I = SSG_2_obj_candidate(1,index_ssg_2_option_I+1);
cons_ssg_2_option_I = SSG_2_cons_candidate(1,index_ssg_2_option_I+1);
w_ssg_2_option_I = SSG_2_w_candidate(:,index_ssg_2_option_I+1);

% print t, the respective s, F(w^(t)) and G(w^(t))
fprintf('Result of SSG_2 with last t in I: iter=%i, s=%f, obj=%f, cons=%f\n', index_ssg_2_option_I, s_ssg_2_option_I, obj_ssg_2_option_I+0, cons_ssg_2_option_I+0);

% Option II: return w^(t) with the minimal value among all F(w^(t))'s with t in I

% select objective values with t in I
SSG_2_obj_candidate_index_set = SSG_2_obj_candidate(SSG_2_index_set+1);
[obj_ssg_2_option_II, index_ssg_2_option_II] = min(SSG_2_obj_candidate_index_set); % get the index in I
SSG_2_index_set_in_cons_candidate = SSG_2_index_set+1; % t in I is 1 less than t in obj and cons candidate
index_ssg_2_option_II = SSG_2_index_set_in_cons_candidate(index_ssg_2_option_II); % get the index in 0,1,...,T
s_ssg_2_option_II = SSG_2_s_candidate(1,index_ssg_2_option_II);
cons_ssg_2_option_II = SSG_2_cons_candidate(1,index_ssg_2_option_II);
w_ssg_2_option_II = SSG_2_w_candidate(:,index_ssg_2_option_II);

% print t, the respective s, F(w^(t)) and G(w^(t))
fprintf('Result of ssg_2 with smallest objective and t in I: iter=%i, s=%f, obj=%f, cons=%f\n', index_ssg_2_option_II-1, s_ssg_2_option_II, obj_ssg_2_option_II+0, cons_ssg_2_option_II+0);

% Option III: return bar(w)^(T) = ( sum_{t in I} eta_t * w^(t) ) / ( sum_{t in I} eta_t )
eta_candidate_index_set = SSG_2_eta_candidate(SSG_2_index_set+1);
w_candidate_index_set = SSG_2_w_candidate(:,SSG_2_index_set+1);
SSG_2_bar_w = w_candidate_index_set*(eta_candidate_index_set') / sum(eta_candidate_index_set);

% compute mean_male = mean( sigmoid ( Ac_male * w - Sc_w ) )
sigmoid_mean_male = mean( 1 ./ ( 1 + exp( Sc_w - Ac(idmalec,:)*SSG_2_bar_w ) ) );
% compute mean_female = mean( sigmoid ( Ac_female * w - Sc_w ) )
sigmoid_mean_female = mean( 1 ./ ( 1 + exp( Sc_w - Ac(idfemalec,:)*SSG_2_bar_w ) ) );
% get the array of objective candidates abs( mean_male - mean_female ) 
% then obtain the largest value of objective F(bar(w)^(T)) with index of s
[obj_ssg_2_option_III, index] = max( abs( sigmoid_mean_male - sigmoid_mean_female ) ); 

% get the value of s
s_ssg_2_option_III = S(1,index); 
% compute the value of constraint G(bar(w)^(T))
cons_ssg_2_option_III = dot( 1 - b.*(A*SSG_2_bar_w), 1 - b.*(A*SSG_2_bar_w) > 0 ) / n - C;
        
fprintf('Result of ssg_2 with weighted w in I: s=%f, obj=%f, cons=%f\n', s_ssg_2_option_III, obj_ssg_2_option_III+0, cons_ssg_2_option_III+0);