function [ RankingLoss,OneError,Coverage,Average_Precision,Macro_AUC,trained_w,trained_b ] = iMLCU( label_train_data,label_train_target,unlabel_train_data,test_data,test_target,C1,C2 )
%iMLCU [2] deals with inductive semi-supervised multi-label learning problem by imposing the max-margin regulatization on unlabeled data.
%
%    Syntax
%
%       [ RankingLoss,OneError,Coverage,Average_Precision,trained_w,trained_b ]=iMLCU( label_train_data,label_train_target,unlabel_train_data,test_data,test_target,C1,C2 )
%
%    Description
%
%       iMLCU takes,
%           label_train_data       - An lxd array, the ith instance of labeled training instance is stored in label_train_data(i,:)
%           label_train_target     - A qxl array, if the ith labeled training instance belongs to the jth class, then label_train_target(j,i) equals +1, otherwise label_train_target(j,i) equals -1
%           unlabel_train_data     - An uxd array, the ith instance of unlabeled training instance is stored in unlabel_train_data(i,:)
%           test_data              - An txd array, the ith instance of testing instance is stored in test_data(i,:)
%           test_target            - A qxt array, if the ith testing instance belongs to the jth class, test_target(j,i) equals +1, otherwise test_target(j,i) equals -1
%           C1                     - The first parameter
%           C2                     - The second parameter
%
%      and returns,
%           RankingLoss      - The ranking loss on testing data as described in [1]
%           OneError         - The one-error on testing data as described in [1]
%           Coverage         - The coverage on testing data as described in [1]
%           Average_Precision- The average precision on testing data as described in [1]
%           trained_w        - An d*q array. The weighed vector of i-th label is stored in trained_w(:,i)
%           trained_b        - An 1*q vector. The bias of i-th label is stored in trained_b(1,i)
%
%  [1] R. E. Schapire, Y. Singer. BoosTexter: A boosting based system for text categorization. Machine Learning, 39(2/3): 135-168, 2000.
%  [2] L. Wu, M.-L. Zhang. Multi-label classification with unlabeled data: An inductive approach. In: Proceedings of the 5th Asian Conference on Machine Learning (ACML'13), Canberra, Australia, 2013.


[num_label_train,dim]=size(label_train_data);
num_unlabel_train = size(unlabel_train_data,1);
num_class = size(label_train_target,1);

for m = 1:num_label_train
    Y(m) = size(label_train_target(label_train_target(:,m) == 1,m),1); % Y(m): number of positive labels of m-th labeled example
    Ybar(m) = num_class - Y(m);
    coef(m,1) = 1/(Y(m)*Ybar(m));
end
coef = C1*coef;   % C1: hyperparameter to be tuned 

%**************************** Get the initial linear models(w and b) on labeled data ********************%
model = cell(num_class,1);   % model for LIBSVM
w = zeros(dim,num_class);
b = zeros(1,num_class);

% initial w and b for each class label on labeled train data
for j = 1:num_class
     model{j,1}  = svmtrain(label_train_target(j,:)',label_train_data,'-t 0 -b 1');
     w(:,j) = model{j,1}.SVs'*model{j,1}.sv_coef; 
     b(1,j) = -model{j,1}.rho;
     if model{j,1}.Label(1) == -1
         w(:,j) = -w(:,j);
         b(1,j) = -b(1,j);
     end
     if size(w(:,j),1) ~= dim
         error('wrong weighted vector w and dimension of w is not same as feature dimension');
     end
end

% store w and b learned from labeled train data
original_w = w;
original_b = b;
%********************************************************************************************************%

Eps_on_labeledTrain = cell(num_label_train,1); % Eps_on_labeledTrain{i,1}: q *q eps matrix of i-th labeled train example
for m = 1:num_label_train
    Eps_on_labeledTrain{m,1} = zeros(num_class);
end

xwb = label_train_data * w;  % l * q matrix
for j = 1:num_class
    xwb(:,j) = xwb(:,j) + b(1,j);
end
unlabel_xwb = unlabel_train_data * w;   % u * q matrix
for j = 1:num_class
    unlabel_xwb(:,j) = unlabel_xwb(:,j) + b(1,j);
end

pre_predictlabel = ones(num_class,num_unlabel_train);  % q * u matrix
opts = optimset('Algorithm','active-set','Display','off','MaxIter',800);
countIteration = 1;   % max iteration number
while 1    
    preValue = valueOfOptimization(w,b,coef,Eps_on_labeledTrain,unlabel_train_data,unlabel_xwb,C2);
    
    for i = 1:num_class   % learn w and b of i-th label
        pre_predictlabel(i,unlabel_xwb(:,i) < 0) = -1;
        
        Eps = [];
        for m = 1:num_label_train
           if label_train_target(i,m) == 1
               epsi = Eps_on_labeledTrain{m,1}(i,label_train_target(:,m) == -1);       
           else
               epsi = Eps_on_labeledTrain{m,1}(i,label_train_target(:,m) == 1);
           end
           Eps = [Eps;epsi']; 
        end
    
        H = zeros((dim+1+size(Eps,1)));
        H(1:dim,1:dim) = eye(dim);
        
        A_C1 = [];  
        A_C2 = [];  
        B_C1 = [];
        F_eps = [];      
        columnFlag = dim+1+1;  
        for m = 1:num_label_train   
            f_eps = [];
            if label_train_target(i,m) == 1
                Am = zeros(Ybar(m),dim+1+size(Eps,1));
                Am(:,1:dim) = repmat(label_train_data(m,:),Ybar(m),1);
                Am(:,dim+1) = 1;
                Am(:,columnFlag:columnFlag+Ybar(m)-1) = eye(Ybar(m));
                Am = (-1)*Am;
            
                Am_eps = zeros(Ybar(m),dim+1+size(Eps,1));
                Am_eps(:,columnFlag:columnFlag+Ybar(m)-1) = eye(Ybar(m));
                Am_eps = (-1)*Am_eps;
            
                Bm = zeros(Ybar(m),1);
                Bm = 1 + xwb(m,label_train_target(:,m) == -1)';
                Bm = Bm * (-1);
            
                f_eps(1:Ybar(m),1) = coef(m,1); 
            
                columnFlag = columnFlag + Ybar(m);
            else
                Am = zeros(Y(m),dim+1+size(Eps,1));
                Am(:,1:dim) = repmat(label_train_data(m,:),Y(m),1);
                Am(:,dim+1) = 1;
                Am(:,columnFlag:columnFlag+Y(m)-1) = eye(Y(m))*(-1);
            
                Am_eps = zeros(Y(m),dim+1+size(Eps,1));
                Am_eps(:,columnFlag:columnFlag+Y(m)-1) = eye(Y(m));
                Am_eps = (-1)*Am_eps;
            
                Bm = zeros(Y(m),1);
                Bm = xwb(m,label_train_target(:,m) == 1)' - 1;
            
                f_eps(1:Y(m),1) = coef(m,1);
            
                columnFlag = columnFlag + Y(m);
            end
            A_C1 = [A_C1;Am];
            A_C2 = [A_C2;Am_eps];
            B_C1 = [B_C1;Bm];
            F_eps = [F_eps;f_eps];
        end
        A = [A_C1;A_C2];
        B_C2 = zeros(size(Eps,1),1);
        B = [B_C1;B_C2];
   
        Aeq = [(1/num_unlabel_train)*sum(unlabel_train_data),1,zeros(1,size(Eps,1))];
        beq = (1/num_label_train)*sum(label_train_target(i,:));
    
        F_wb = zeros(dim+1,1);
        for n = 1:num_unlabel_train
            if pre_predictlabel(i,n) == 1
                if unlabel_xwb(n,i) <= -1
                    f_wb = (-2)*[unlabel_train_data(n,:),1]'*C2; 
                else if unlabel_xwb(n,i) < 1
                        f_wb = ((-1)*[unlabel_train_data(n,:),1]')*C2;
                    else
                        f_wb = zeros(dim+1,1);
                    end
                end
            else
                if unlabel_xwb(n,i) >= 1
                    f_wb = 2*[unlabel_train_data(n,:),1]'*C2;
                else if unlabel_xwb(n,i) > -1
                      f_wb = ([unlabel_train_data(n,:),1]')*C2;
                    else
                        f_wb = zeros(dim+1,1);
                    end
                end
            end
            F_wb = F_wb + f_wb;
        end
        F = [F_wb;F_eps];
        
        x0 = [w(:,i);b(1,i);Eps];  % x0: (dim+1+size(Eps)) * 1 vector
        %[X,fval,exitflag,output,lambda] = quadprog(H,F,A,B,Aeq,beq,[],[],x0,opts);
        X = quadprog(H,F,A,B,Aeq,beq,[],[],x0,opts);
        w(:,i) = X(1:dim);
        b(1,i) = X(dim+1);
        
        Eps = X(dim+2:end);
        unlabel_xwb(:,i) = unlabel_train_data * w(:,i) + b(1,i);  % update the outputs of unlabeled data on i-th w and b
        xwb(:,i) = label_train_data * w(:,i) + b(1,i);
        
        colIndex = 1;
        for m = 1:num_label_train   % update Eps_on_labeledTrain
           if label_train_target(i,m) == 1
               Eps_on_labeledTrain{m,1}(i,label_train_target(:,m) == -1) = Eps(colIndex:colIndex+Ybar(m)-1);
               colIndex = colIndex + Ybar(m);
           else
               Eps_on_labeledTrain{m,1}(i,label_train_target(:,m) == 1) = Eps(colIndex:colIndex+Y(m)-1);
               colIndex = colIndex + Y(m);
           end
           Eps_on_labeledTrain{m,1}(:,i) = Eps_on_labeledTrain{m,1}(i,:);  % Eps_on_labeledTrain is symmetric matrix
        end        
        
    end
    nextValue = valueOfOptimization(w,b,coef,Eps_on_labeledTrain,unlabel_train_data,unlabel_xwb,C2);
    wb_on_everyIter{countIteration,1} = w;
    wb_on_everyIter{countIteration,2} = b;
    Value(1,countIteration) = preValue;
    Value(1,countIteration+1) = nextValue;
    
    if preValue > nextValue && (preValue - nextValue)/nextValue < 0.01
        disp('optimization converge at ratio difference 0.01');
        break;   % break while loop      
    end
    if countIteration < 10
        countIteration = countIteration + 1;
        continue;
    else
        disp('exceed max iteration');
        break;   % break while loop
    end
    
end  % end while loop
disp(['countIteration: ',num2str(countIteration)]);

% choose the w and b which generate the minimum value of the optimization formula
[~,index] = sort(Value,'ascend');   
if index(1) == 1
    trained_w = original_w;
    trained_b = original_b;
else
    trained_w = wb_on_everyIter{index(1)-1,1};
    trained_b = wb_on_everyIter{index(1)-1,2};
end

test_xwb = test_data * trained_w;
for i = 1:num_class
    test_xwb(:,i) = test_xwb(:,i) + trained_b(1,i);
end
Pre_Label = ones(size(test_xwb));
Pre_Label(test_xwb < 0) = -1;
Pre_Label = Pre_Label';
Outputs = test_xwb';

% compute the metrics on the test dataset
RankingLoss = Ranking_loss(Outputs,test_target);
OneError = One_error(Outputs,test_target);
Coverage = coverage(Outputs,test_target);
Average_Precision = Average_precision(Outputs,test_target);
Macro_AUC = macro_auc(Outputs,test_target);

end  % end function iMLCU

