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

% prepare an 1*(numiterssg_1_ub+1) zero array to store the values of objective in each iteration
SSG_1_obj_candidate = zeros(1, numiterssg_1_ub+1);
% prepare an 1*(numiterssg_1_ub+1) zero array to store the values of s in each iteration
SSG_1_s_candidate = zeros(1, numiterssg_1_ub+1);
% prepare an 1*(numiterssg_1_ub+1) zero array to store the values of constraint in each iteration
SSG_1_cons_candidate = zeros(1, numiterssg_1_ub+1);
% prepare a d*(numiterssg_1_ub+1) zero matrix to store the iterates in each iteration by column
SSG_1_w_candidate = zeros(d, numiterssg_1_ub+1);

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

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

% initialize I in SSG method
SSG_1_index_set = [];
% initialize CPU time counting in SSG method
SSG_1_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_1_w_candidate(:,1) = w_init; 
SSG_1_time_set(length(SSG_1_time_set)) = SSG_1_time_set(length(SSG_1_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_1_s_candidate(1,1) = S(1,index); 
% get the value of objective F(w^(0))
SSG_1_obj_candidate(1,1) = obj; 

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

fprintf('Iter=%i, s=%f, obj=%f, cons=%f\n', 0, SSG_1_s_candidate(1,1), SSG_1_obj_candidate(1,1)+0, SSG_1_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 ) / delta;
% epsilon = 0.001;

for t = 1:numiterssg_1_ub
   
    tic; % for recording CPU runtime

    % epsilon_t = epsilon^2 * (rho_hat - rho) / (1 + Lambda); 
    ssg_1_epsilon_t = 1e-5; % select from {1e-6, 2e-6, 5e-6, 1e-5}

    if SSG_1_cons_candidate(1,t)+0 < ssg_1_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_1_obj_subgrad_candidate(:,t) = obj_subgrad;
        
        % eta_t = epsilon^2 * (rho_hat - rho) / ( (1 + Lambda) * M^2 ); % compute eta_t
        ssg_1_eta_t = 1e-3; % select from {2e-4, 5e-4, 1e-3, 2e-3}
        SSG_1_eta_candidate(1,t) = ssg_1_eta_t;
        
        w = w - ssg_1_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_1_w_candidate(:,t+1) = w;

        SSG_1_index_set = [SSG_1_index_set, t-1]; % add t to I
        SSG_1_time_set = [SSG_1_time_set, SSG_1_time_set(length(SSG_1_time_set))+toc]; % record CPU runtime

        step_count_ssg_1 = step_count_ssg_1 + 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_1_cons_subgrad_candidate(:,t) = cons_subgrad;
        
        % eta_t = epsilon^2 * (rho_hat - rho) / ( (1 + Lambda) * M^2 ); % compute eta_t
        ssg_1_eta_t = 1e-3; % select from {2e-4, 5e-4, 1e-3, 2e-3}
        SSG_1_eta_candidate(1,t) = ssg_1_eta_t;
        
        w = w - ssg_1_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_1_w_candidate(:,t+1) = w;

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

        step_count_ssg_1 = step_count_ssg_1 + 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_1_s_candidate(1,t+1) = S(1,index); 
    % get the value of objective F(w^(t))
    SSG_1_obj_candidate(1,t+1) = obj; 

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

    fprintf('Iter=%i, s=%f, obj=%f, cons=%f\n', t, SSG_1_s_candidate(1,t+1), SSG_1_obj_candidate(1,t+1)+0, SSG_1_cons_candidate(1,t+1)+0);
   
end

% Option I: return w^(t) with the largest t in I
index_ssg_1_option_I = SSG_1_index_set(1,length(SSG_1_index_set));
s_ssg_1_option_I = SSG_1_s_candidate(1,index_ssg_1_option_I+1);
obj_ssg_1_option_I = SSG_1_obj_candidate(1,index_ssg_1_option_I+1);
cons_ssg_1_option_I = SSG_1_cons_candidate(1,index_ssg_1_option_I+1);
w_ssg_1_option_I = SSG_1_w_candidate(:,index_ssg_1_option_I+1);

% print t, the respective s, F(w^(t)) and G(w^(t))
fprintf('Result of SSG_1 with last t in I: iter=%i, s=%f, obj=%f, cons=%f\n', index_ssg_1_option_I, s_ssg_1_option_I, obj_ssg_1_option_I+0, cons_ssg_1_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_1_obj_candidate_index_set = SSG_1_obj_candidate(SSG_1_index_set+1);
[obj_ssg_1_option_II, index_ssg_1_option_II] = min(SSG_1_obj_candidate_index_set); % get the index in I
SSG_1_index_set_in_cons_candidate = SSG_1_index_set+1; % t in I is 1 less than t in obj and cons candidate
index_ssg_1_option_II = SSG_1_index_set_in_cons_candidate(index_ssg_1_option_II); % get the index in 0,1,...,T
s_ssg_1_option_II = SSG_1_s_candidate(1,index_ssg_1_option_II);
cons_ssg_1_option_II = SSG_1_cons_candidate(1,index_ssg_1_option_II);
w_ssg_1_option_II = SSG_1_w_candidate(:,index_ssg_1_option_II);

% print t, the respective s, F(w^(t)) and G(w^(t))
fprintf('Result of SSG_1 with smallest objective and t in I: iter=%i, s=%f, obj=%f, cons=%f\n', index_ssg_1_option_II-1, s_ssg_1_option_II, obj_ssg_1_option_II+0, cons_ssg_1_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_1_eta_candidate(SSG_1_index_set+1);
w_candidate_index_set = SSG_1_w_candidate(:,SSG_1_index_set+1);
ssg_1_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_1_bar_w ) ) );
% compute mean_female = mean( sigmoid ( Ac_female * w - Sc_w ) )
sigmoid_mean_female = mean( 1 ./ ( 1 + exp( Sc_w - Ac(idfemalec,:)*ssg_1_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_1_option_III, index] = max( abs( sigmoid_mean_male - sigmoid_mean_female ) ); 

% get the value of s
s_ssg_1_option_III = S(1,index); 
% compute the value of constraint G(bar(w)^(T))
cons_ssg_1_option_III = dot( 1 - b.*(A*ssg_1_bar_w), 1 - b.*(A*ssg_1_bar_w) > 0 ) / n - C;

fprintf('Result of SSG_1 with weighted w in I: s=%f, obj=%f, cons=%f\n', s_ssg_1_option_III, obj_ssg_1_option_III+0, cons_ssg_1_option_III+0);