function measure = deterministic_weakly_convex_stationarity_measure(A, b, Ac, idmalec, idfemalec, nummalec, numfemalec, w, rho_hat, SCAD_coefficient, kappa)
    
    [n,~] = size(A);

    % prepare an upper bound of number of iterations for the algorithm
    numiter_ub = 2500;
    % initialize the counting on the actual number of iterations in the algorithm
    step_count = 0;
    % initialize the tolerance for the stopping criteria
    tol = 0;
    % prepare an indicator to determine whether break the iteration or not
    break_indicator = 0;

    % w_0 is inherited from the current iteration, which is d*1 double
    w_0 = w;
    w_old = w; 
    
    % initialize I in the algorithm
    index_set = [];
    
    gamma = 1e-5; 

    for t = 1:numiter_ub
   
        epsilon_t = 1e-5; 

        % 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)
        cons = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa + rho_hat / 2 * norm(w-w_0, 2)^2 - gamma;

        if cons+0 < 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) ...
                          + rho_hat * ( w - w_0 ); % key difference

            eta_t = 7.5e-4; 
            
            w_old = w;
            w = w - 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

            step_count = step_count + 1;
            index_set = [index_set, t-1]; % add t to I

            % stopping criteria
            if length(index_set) > 1 % start to compare if there are two elements in I

                % compare two adjacent iterates in I
                if norm(w_old-w, 2) <= tol
                    break_indicator = 1; % change the indicator
                end
            end

        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) ...
                           + rho_hat * ( w - w_0 ); % key difference

            eta_t = 7.5e-4; 

            w = w - 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

            step_count = step_count + 1;

        end

        % determine whether break or not according to the indicator
        if break_indicator == 1
            break
        end
        
        % compute the value of objective F(w^(t+1))
        obj = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n + SCAD_coefficient * SCAD(w) ...
              + rho_hat / 2 * norm(w-w_0, 2)^2; 
        
        % 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))
        cons = abs( sigmoid_mean_male - sigmoid_mean_female ) - kappa ...
               + rho_hat / 2 * norm(w-w_0, 2)^2 - gamma;
        
        fprintf('Iter=%i, obj=%f, cons=%f\n', t, obj+0, cons+0);
    end

    measure = norm(w-w_0, 2);
    fprintf('compute the stationarity with %i steps to get the stationarity=%f\n', step_count, measure);
end