% prepare an upper bound of number of iterations for the algorithm
numiterippconex_ub = 600;
% prepare an upper bound of number of sub-iterations for the algorithm
numitersub_ub = 600;
% initialize the counting on the actual number of iterations in the algorithm
step_count_ippconex = 0;

% prepare a d*(numiterippconex_ub+1) zero matrix to store the iterates in each iteration by column
IPP_ConEx_w_candidate = zeros(d, numiterippconex_ub+1);
% prepare an 1*(numiterippconex_ub+1) zero array to store the values of objective in each iteration
IPP_ConEx_obj_candidate = zeros(1, numiterippconex_ub+1);
% prepare an 1*(numiterippconex_ub+1) zero array to store the values of constraint in each iteration
IPP_ConEx_cons_candidate = zeros(1, numiterippconex_ub+1);

% initialize CPU time counting in IPP-ConEx method
IPP_ConEx_time_set = [0];

% w_init is inherited from the initialization, which is d*1 double
w = w_init;
IPP_ConEx_w_candidate(:,1) = w_init;

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

% 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 objective G(w^(0))
IPP_ConEx_cons_candidate(1,1) = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa;

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

for t = 1:numiterippconex_ub
 
    % prepare an 1*(numitersub_ub+1) zero array to store the values of objective in each subiteration
    obj_candidate = zeros(1, numitersub_ub+1);
    % prepare an 1*(numitersub_ub+1) zero array to store the values of constraints in each subiteration
    cons_candidate = zeros(1, numitersub_ub+1);
    % prepare an 1*(numitersub_ub+1) zero array to store the values of linear approximation 
    % of constraint in each subiteration
    l_cons_candidate = zeros(1, numitersub_ub+1);
    % prepare a d*(numitersub_ub+1) zero matrix to store the iterates in each subiteration by column
    w_candidate = zeros(d, numitersub_ub+1);
    % prepare an 1*(numitersub_ub+1) zero array to store the dual iterates in each subiteration
    z_candidate = zeros(1, numitersub_ub+1);
    
    % prepare a d*(numitersub_ub) zero matrix to store the subgradients of constraint in each iteration by column
    cons_subgrad_candidate = zeros(d, numitersub_ub);

    % initialize I in subiteration of ippconex method
    index_set = [];

    tic; % for recording CPU runtime
    % w^(t,0) is from the initialization or the previous iteration
    w_candidate(:,1) = w; 
    % z^(t,0) = 0
    z = 0;
    z_candidate(1,1) = 0; 
    IPP_ConEx_time_set(length(IPP_ConEx_time_set)) = IPP_ConEx_time_set(length(IPP_ConEx_time_set)) + toc; % record CPU runtime

    % compute the value of objective
    % F_t(w) = F(w) + rho_hat / 2 * norm(w-w^(t,0), 2)^2 on w = w^(t,0)
    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_t(w) = G(w) + rho_hat / 2 * norm(w-w^(t,0), 2)^2 on w = w^(t,0)
    cons_candidate(1,1) = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa;
    % get the value of linear apprximation l_G(w^(t,0)) = G_t(w^(t,0))
    l_cons_candidate(1,1) = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa;
    IPP_ConEx_time_set(length(IPP_ConEx_time_set)) = IPP_ConEx_time_set(length(IPP_ConEx_time_set)) + toc; % record CPU runtime

    fprintf('Iter=%i, Subiter=%i, obj=%f, cons=%f\n', t-1, 0, obj_candidate(1,1)+0, cons_candidate(1,1)+0);
    % add +0 to avoid the error: using fprintf. Function is not defined for sparse inputs
    
    substep_count_ippconex = 0;

    for k = 1:numitersub_ub
        
        tic; % for recording CPU runtime

        if k == 1
            s_k = l_cons_candidate(1,k);
        else
            l_cons_candidate(1,k) = cons_candidate(1,k-1) + dot( cons_subgrad_candidate(:,k-1), w - w_candidate(:,k-1) );
            s_k = ( 2 - 1 / k ) * l_cons_candidate(1,k) - k / (k+1) * l_cons_candidate(1,k-1);
        end
        
        % set tau_k
        ippconex_tau_k = 500 * k; % select from {50, 100, 200, 500}
        z = max(z + ippconex_tau_k*s_k, 0); % update z^(k)
        z_candidate(1,k+1) = z; 

        % compute the subgradients for both objective and constraint for the update
        % 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) ...
                      + rho_hat * ( w - w_candidate(:,1) );

        % compute constraint subgradient by calling the function "constraint_subgradient"
        cons_subgrad = constraint_subgradient(Ac, idmalec, idfemalec, nummalec, numfemalec, w) ...
                       + rho_hat * ( w - w_candidate(:,1) );
        cons_subgrad_candidate(:,k) = cons_subgrad;

        % set eta_k
        ippconex_eta_k = 5e-3 / k; % select from {5e-3, 1e-2, 2e-2, 5e-2}
        w = w - ippconex_eta_k * ( obj_subgrad + z * cons_subgrad); % update w^(k)
        % w = min( 0.5*D_X / norm(w,2), 1 ) * w; % projection to the ball with radius as 0.5*D_X
        w_candidate(:,k+1) = w;
        
        IPP_ConEx_time_set = [IPP_ConEx_time_set, IPP_ConEx_time_set(length(IPP_ConEx_time_set))+toc]; % record CPU runtime
        
        step_count_ippconex = step_count_ippconex + 1;
        substep_count_ippconex = substep_count_ippconex + 1;
    
        % compute the value of objective
        % F_t(w^(k)) = F(w^(t,k)) + rho_hat / 2 * norm(w^(t,k)-w^(t,0), 2)^2
        obj_candidate(1,k+1) = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n + SCAD_coefficient * SCAD(w) ...
                               + rho_hat / 2 * norm(w-w_candidate(:,1), 2)^2;
        
        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_t(w^(k)) = G(w^(t,k)) + rho_hat / 2 * norm(w^(t,k)-w^(t,0), 2)^2
        cons_candidate(1,k+1) = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa ...
                                + rho_hat / 2 * norm(w-w_candidate(:,1), 2)^2;
        IPP_ConEx_time_set(length(IPP_ConEx_time_set)) = IPP_ConEx_time_set(length(IPP_ConEx_time_set)) + toc; % record CPU runtime

        fprintf('Iter=%i, Subiter=%i, obj=%f, cons=%f\n', t-1, k, obj_candidate(1,k+1)+0, cons_candidate(1,k+1)+0);
     
    end
    
    tic; % for recording CPU runtime
    % return bar(w)^(K) = ( sum_{k} (k+1) * w^(k) ) / ( sum_{k} (k+1) )
    k_candidate = 1:(substep_count_ippconex+1);
    bar_w = w_candidate*(k_candidate') / sum(k_candidate);
    
    % w = w^(t) = w^(t,0) = bar(w)^(K)
    w = bar_w;
    IPP_ConEx_w_candidate(:,t+1) = w; 
    IPP_ConEx_time_set(length(IPP_ConEx_time_set)) = IPP_ConEx_time_set(length(IPP_ConEx_time_set)) + toc; % record CPU runtime

    % compute the value of constraint F(w^(t))
    obj = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n + SCAD_coefficient * SCAD(w); 
    IPP_ConEx_obj_candidate(1,t+1) = obj;
    
    % 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 objective G(w^(t))
    cons = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa;
    IPP_ConEx_cons_candidate(1,t+1) = cons;

    fprintf('Result of Iter=%i in ippconex: obj=%f, cons=%f\n', t-1, obj+0, cons+0);

end

% return w^(T) = bar(w)^(K) in the last iteration
fprintf('Result of ippconex with last iteration: iter=%i, obj=%f, cons=%f\n', step_count_ippconex, obj+0, cons+0);