% Part I: prepare the dataset

load('a9a.mat'); % a9a.mat contains sparse double matrix A with size 32561*123, 
                 % sparse double matrix Ac with size 16281*123 
                 % and double vector b with size 32561*1.

% On the page https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html#a9a, 
% download a9a as a9a.txt and a9a.t as a9a.t.txt.
% load('a9a.mat');
% [b_1,  A_1] = libsvmread('C:\Program Files\MATLAB\R2021a\toolbox\libsvm-3.3\a9a.txt');
% [bc_1,  Ac_1] = libsvmread('C:\Program Files\MATLAB\R2021a\toolbox\libsvm-3.3\a9a.t.txt');
% isequal(A,  A_1)
% ans = logical  1
% isequal(b, b_1)
% ans = logical  1
% isequal( full(Ac), [ bc_1, full(Ac_1) ] )
% ans = logical  1

% i) [A, b] is for constructing the constaint.

[n,d] = size(A); % number of instances n = 32561, feature dimension d = 123
numclass = length(unique(b)); % b is vector of class labels.
% There are two labels: 1 and -1, so numclass = 2.
% To access an element in the sparse double matrix A, use full(A(i,j)) to 
% get the value of (i,j)-th element.

sv = full(A(:,72)); % sv = 0 for "male" and = 1 for "female"
idmale = find(sv==0); % indexes in [n] of instances for "male"
idfemale = find(sv==1); % indexes in [n] of instances for "female"

nummale = sum(sv==0); % number of instances for "male" = 21790
numfemale = sum(sv==1); % number of instances for "female" = 10771

numpos = sum(b==1); % number of instances with positive, i.e. label 1 = 7841
numneg = sum(b==-1); % number of instances with negative, i.e. label -1 = 24720

numposmale = sum(b==1 & sv==0); % number of instances for "male" and with label 1 = 6662
numnegmale = sum(b==-1 & sv==0); % number of instances for "male" and with label -1 = 15128
numposfemale = sum(b==1 & sv==1); % number of instances for "female" and with label 1 = 1179
numnegfemale = sum(b==-1 & sv==1); % number of instances for "female" and with label -1 = 9592

% ii) [Ac, bc] is for constructing the objective.

bc = full(Ac(:,1));
Ac = Ac(:,2:size(Ac, 2));
Ac = [Ac,repelem(0,size(Ac, 1),1)];

[nc,dc] = size(Ac); % number of instances n = 16281, feature dimension d = 123

numclassc = length(unique(bc)); % bc is vector of class labels.
% There are two labels: 1 and -1, so numclass = 2.

% sensitive variable, abbreviated by sv, concerns the fairness.
svc = full(Ac(:,72)); % sv = 0 for "male" and = 1 for "female"
idmalec = find(svc==0); % indexes in [nc] of instances for "male"
idfemalec = find(svc==1); % indexes in [nc] of instances for "female"

nummalec = sum(svc==0); % number of instances for "male" = 10860
numfemalec = sum(svc==1); % number of instances for "female" = 5421

numposc = sum(bc==1); % number of instances with positive, i.e. label 1 = 3846
numnegc = sum(bc==-1); % number of instances with negative, i.e. label -1 = 12435

numposmalec = sum(bc==1 & svc==0); % number of instances for "male" and with label 1 = 3256
numnegmalec = sum(bc==-1 & svc==0); % number of instances for "male" and with label -1 = 7604
numposfemalec = sum(bc==1 & svc==1); % number of instances for "female" and with label 1 = 590
numnegfemalec = sum(bc==-1 & svc==1); % number of instances for "female" and with label -1 = 4831

% Part II: initialize w^(0) by zero vector
beta = ( mean(vecnorm(full(Ac(idmalec,:)),2,2).^2) + mean(vecnorm(full(Ac(idfemalec,:)),2,2).^2) ) / 4;
SCAD_coefficient = 0.2;
rho = max(beta, 2*SCAD_coefficient);
rho_hat = max(rho, 1) * 1; % select from {1, 1.5, 2} 

M = max( sqrt( eigs( ((full(A)').*b') * (b.*full(A)), 1 ) ) / n, mean(vecnorm(full(Ac),2,2)) );

kappa = 0.005; % 0.2

w_init = zeros(d,1); 
% D_X = 10;

% Part III: solve the constrained problem by three methods

% apply the SSG method by constant step-size
% to solve the constrained problem
SSG_deterministic_weakly_convex_minimization
w_SSG_result = w;

% apply the IPP-SSG method to solve the constrained problem
IPP_SSG_deterministic_weakly_convex_minimization 
w_IPP_SSG_result = w;
 
% apply the IPP-ConEx method to solve the constrained problem
IPP_ConEx_deterministic_weakly_convex_minimization 
w_IPP_ConEx_result = w;

% Part IV: visualize the results

title_name = "a9a";
% deterministic_weakly_convex_visualization