function measure = deterministic_convex_stationarity_measure_batch(A, b, Ac, idmalec, idfemalec, nummalec, numfemalec, w, rho_hat, C, D_X)
    
    [n,~] = size(A);
    [~,num_measures] = size(w); % each COLUMN of w is a solution to compute stationarity

    % 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 = zeros(num_measures,1);

    % 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 = cell(num_measures, 1);
    index_set_length = zeros(num_measures, 1);

    % 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
    % Compute Sc_w for each w, Sc_w is a matrix, each column correspond a w
    Sc_w = ones(size(S,2),1)* min(Ac*w) + S'*( max(Ac*w) - min(Ac*w) ); 
    
    for t = 1:numiter_ub
   
        epsilon_t = 1e-5; 
        cons = dot( 1 - b.*(A*w), 1 - b.*(A*w) > 0 ) / n - C;
        % generate id of w that needs to go along subgradient of objective
        wid_obj = find(cons+0 < epsilon_t);
        wid_obj = setdiff(wid_obj, find(break_indicator));
        % generate id of w that needs to go along subgradient of constraint
        wid_con = find(cons+0 >= epsilon_t);
        wid_con = setdiff(wid_con, find(break_indicator));
        
        % compute objevtive subgradient by calling the function "objective_subgradient_batch"
        if length(wid_obj) > 0
            obj_subgrad = objective_subgradient_batch(Ac, idmalec, idfemalec, nummalec, numfemalec, w(:, wid_obj), Sc_w(:, wid_obj)) ...
                          + rho_hat * ( w(:, wid_obj) - w_0(:, wid_obj) ); % key difference

            eta_t = 2e-3; 
            w_old(:, wid_obj) = w(:, wid_obj);
            w(:, wid_obj) = w(:, wid_obj) - eta_t * obj_subgrad; % update w^(t)
            
            for i = wid_obj
                index_set{i} = [index_set{i}, t-1]; % add t to I
            end

            % stopping criteria
            if min(index_set_length) > 1 % start to compare if there are two elements in I
                break_indicator( wid_obj( vecnorm(w_old(:, wid_obj)-w(:, wid_obj), 2) <=tol ) ) = 1;
            end
        end        

        if length(wid_con) > 0
            % compute the constraint subgradient as the same as that in hingeloss_minimization
            cons_subgrad = - ( (A').*b' ) * ( 1 - b.*(A*w(:, wid_con)) > 0 ) / n;
            
            eta_t = 2e-3; 
            w(:, wid_con) = w(:, wid_con) - eta_t * cons_subgrad; % update w^(t)
        end                  

        w = min( 0.5*D_X ./ vecnorm(w,2), 1 ) .* w; % projection to the ball with radius as 0.5*D_X

        % determine whether break or not according to the indicator
        if min(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);
        fprintf('Iter=%i\n', t);
    end

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