function measure = deterministic_weakly_convex_stationarity_measure_batch(A, b, Ac, idmalec, idfemalec, nummalec, numfemalec, w, rho_hat, SCAD_coefficient, kappa)
    
    [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);
    
    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 * vecnorm(w-w_0, 2).^2 - gamma;
        
        % 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_cons = find(cons+0 >= epsilon_t);
        wid_cons = setdiff(wid_cons, find(break_indicator));

        if length(wid_obj) > 0 % perform subgradient on objective function for wid_obj in w
            
            % compute the objective subgradient as a sum of that of hinge_loss and that of SCAD
            obj_subgrad = - ( (A').*b' ) * ( 1 - b.*(A*w(:, wid_obj)) > 0 ) / n ...
                          + SCAD_coefficient * SCAD_subgradient(w(:, wid_obj)) ...
                          + rho_hat * ( w(:, wid_obj) - w_0(:, wid_obj) ); % key difference

            eta_t = 2e-4; 
            
            w_old(:, wid_obj) = w(:, wid_obj);
            w(:, wid_obj) = w(:, wid_obj) - eta_t * obj_subgrad; % update w^(t)
            % w = min( 0.5*D_X ./ vecnorm(w,2), 1 ) .* w; % projection to the ball with radius as 0.5*D_X

            step_count = step_count + 1;
            for i = wid_obj
                index_set{i} = [index_set{i}, t-1]; % add t to I
            end
        end

        if length(wid_cons) > 0 % perform subgradient on constraint function for wid_cons in w
       
            % compute constraint subgradient by calling the function "constraint_subgradient"
            cons_subgrad = constraint_subgradient_batch(Ac, idmalec, idfemalec, nummalec, numfemalec, w(:, wid_cons)) ...
                           + rho_hat * ( w(:, wid_cons) - w_0(:, wid_cons) ); % key difference

            eta_t = 2e-4; 

            w(:, wid_cons) = w(:, wid_cons) - eta_t * cons_subgrad; % update w^(t)
            % w = min( 0.5*D_X ./ vecnorm(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
        
        fprintf('Iter=%i\n', t);
    end

    measure = vecnorm(w-w_0, 2);
end