function measure = deterministic_convex_stationarity_measure(A, b, Ac, idmalec, idfemalec, nummalec, numfemalec, w, rho_hat, C, D_X)
    
    [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 = [];

    % 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; 
    
    for t = 1:numiter_ub
   
        epsilon_t = 1e-5; 
        cons = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n - C;
        if cons+0 < 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) ...
                          + rho_hat * ( w - w_0 ); % key difference

            eta_t = 2e-3; 
            
            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 the constraint subgradient as the same as that in hingeloss_minimization
            cons_subgrad = - ( (A').*b' ) * ( 1 - b.*(A*w) > 0 ) / n;

            eta_t = 2e-3; 

            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 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, ~] = max( abs( sigmoid_mean_male - sigmoid_mean_female ) ); 
        

        % obj = obj + rho_hat / 2 * norm(w-w_0, 2)^2;
        % cons = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n - C;
        % 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