% prepare an upper bound of number of iterations for the algorithm
numiterssg_ub = 360000;
% initialize the counting on the actual number of iterations in the algorithm
step_count_ssg = 0;

% prepare a d*(numiterssg_ub+1) zero matrix to store the iterates in each iteration by column
SSG_w_candidate = zeros(d, numiterssg_ub+1);
% prepare an 1*(numiterssg_ub+1) zero array to store the values of objective in each iteration
SSG_obj_candidate = zeros(1, numiterssg_ub+1);
% prepare an 1*(numiterssg_ub+1) zero array to store the values of constraint in each iteration
SSG_cons_candidate = zeros(1, numiterssg_ub+1);
% prepare an 1*(numiterssg_ub) zero array to store the step-sizes in each iteration
SSG_eta_candidate = zeros(1, numiterssg_ub);

% initialize I in ssg method
SSG_index_set = [];
% initialize CPU time counting in SSG method
SSG_time_set = [0];

% w_init is inherited from the initialization, which is d*1 double
tic; % for recording CPU runtime
w = w_init;
SSG_w_candidate(:,1) = w_init; 
SSG_time_set(length(SSG_time_set)) = SSG_time_set(length(SSG_time_set)) + toc; % record CPU runtime

% compute the value of objective F(w^(0))
SSG_obj_candidate(1,1) = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n + SCAD_coefficient * SCAD(w); 

tic; % for recording CPU runtime
% compute mean_male = mean( sigmoid ( Ac_male * w ) )
sigmoid_mean_male = mean( 1 ./ ( 1 + exp( - Ac(idmalec,:)*w ) ) );
% compute mean_female = mean( sigmoid ( Ac_female * w ) )
sigmoid_mean_female = mean( 1 ./ ( 1 + exp( - Ac(idfemalec,:)*w ) ) );
% get the value of constraint G(w^(0))
SSG_cons_candidate(1,1) = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa;
SSG_time_set(length(SSG_time_set)) = SSG_time_set(length(SSG_time_set)) + toc; % record CPU runtime

fprintf('Iter=%i, obj=%f, cons=%f\n', 0, SSG_obj_candidate(1,1)+0, SSG_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
c = 1;

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

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

    if SSG_cons_candidate(1,t)+0 < ssg_epsilon_t % perform subgradient on objective function

        % compute the objective subgradient as a sum of that of hinge_loss and that of SCAD
        obj_subgrad = - ( (A').*b' ) * ( 1 - b.*(A*w) > 0 ) / n + SCAD_coefficient * SCAD_subgradient(w);
        
        % eta_t = epsilon^2 * (rho_hat - rho) / ( (1 + Lambda) * M^2 ); % compute eta_t
        ssg_eta_t = 5e-4; % select from {1e-4, 2e-4, 5e-4, 7.5e-4}
        SSG_eta_candidate(1,t) = ssg_eta_t;
        
        w = w - ssg_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_w_candidate(:,t+1) = w;
        
        SSG_index_set = [SSG_index_set, t-1]; % add t to I
        SSG_time_set = [SSG_time_set, SSG_time_set(length(SSG_time_set))+toc]; % record CPU runtime

        step_count_ssg = step_count_ssg + 1;

    else % perform subgradient on constraint function

        % compute constraint subgradient by calling the function "constraint_subgradient"
        cons_subgrad = constraint_subgradient(Ac, idmalec, idfemalec, nummalec, numfemalec, w);
        
        % eta_t = epsilon^2 * (rho_hat - rho) / ( (1 + Lambda) * M^2 ); % compute eta_t
        ssg_eta_t = c * SSG_cons_candidate(1,t) / ( norm(cons_subgrad,2)^2 ); % tune here
        SSG_eta_candidate(1,t) = ssg_eta_t;
        
        w = w - ssg_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_w_candidate(:,t+1) = w;
        
        SSG_time_set = [SSG_time_set, SSG_time_set(length(SSG_time_set))+toc]; % record CPU runtime

        step_count_ssg = step_count_ssg + 1;
        
    end

    % compute the value of objective F(w^(t+1))
    SSG_obj_candidate(1,t+1) = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n + SCAD_coefficient * SCAD(w); 
    
    tic; % for recording CPU runtime
    % compute mean_male = mean( sigmoid ( Ac_male * w ) )
    sigmoid_mean_male = mean( 1 ./ ( 1 + exp( - Ac(idmalec,:)*w ) ) );
    % compute mean_female = mean( sigmoid ( Ac_female * w ) )
    sigmoid_mean_female = mean( 1 ./ ( 1 + exp( - Ac(idfemalec,:)*w ) ) );
    % get the value of constraint G(w^(t+1))
    SSG_cons_candidate(1,t+1) = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa;
    SSG_time_set = [SSG_time_set, SSG_time_set(length(SSG_time_set))+toc]; % record CPU runtime

    fprintf('Iter=%i, obj=%f, cons=%f\n', t, SSG_obj_candidate(1,t+1)+0, SSG_cons_candidate(1,t+1)+0);
   
end