function [M]=COMMU(train_data,train_target,alpha,C,threshold)
%% This function is COMMU method
% INPUTS:
% train_data   (num_data, num_dim) training data
% train_target (num_data, num_label) labels of training data
% alpha  the balancing parameter in Eq.3
% C  the cost parameter in objective function (Eq.6)
% threshold the thresholding parameter in Eq.4
% OUTPUTS:
% M (dim, dim) the learned metric matrix
% For each marker, > 0 is positive, = 0 is negative, Pcenter is normal mean, and Ncenter is negative mean. 

[num_data,num_dim] = size(train_data);
num_label = size(train_target,2);
Pcenter = zeros(num_dim,num_label);
Ncenter = zeros(num_dim,num_label);
for j = 1:num_label
    Pindex = find(train_target(:,j) > 0);
    Nindex = find(train_target(:,j) <= 0);
    Pcenter(:,j) = mean(train_data(Pindex,:))';
    Ncenter(:,j) = mean(train_data(Nindex,:))';
end
B = Pcenter - Ncenter;
B(isnan(B)) = 0;
Bcorr = corr(B);
Bcorr(isnan(Bcorr)) = 0;
Bcorr = abs(Bcorr);
B = mapminmax(B',0,1)';% normalization 

%% The graph structure of label space is constructed. 
% The co-occurrence frequency of markers is taken as the edge, 
% and the adjacent matrix A is symmetrical. It only needs to calculate the first half. 
A = zeros(num_label,num_label);
Anorm = zeros(num_label,num_label);
for j = 1:num_label
    for m = j:num_label
        A(j,m) = size(find(train_target(:,j) == 1 & train_target(:,m) == 1),1); 
        if j ~= m
            A(m,j) = A(j,m);
        end
    end
   
end

Lpnum = diag(A); % number of positive instance for the L-th label

for j = 1:num_label     
    if(Lpnum(j) > 0)
       Anorm(j,:) = A(j,:) / Lpnum(j);
    else
       Anorm(j,:) = A(j,:);
    end
end

% initial parameters:  w = 1/q
w = 1 / num_label * ones(1,num_label);
beta = 1 - alpha;

% Similar constraint set: S; 
% Dissimilar constraint set: R
num_similar = 3;
num_dissimilar = 3;
G = alpha * Anorm + beta * Bcorr;% example_based similarity

Sneighbor = zeros(num_data,num_similar);
Rneighbor = zeros(num_data,num_dissimilar);

valid_index = zeros(1,num_data-1);
select_matrix = ones(1,num_data-1);
for i = 1:num_data
    pp = 1;
    for j = 1:num_data
        if i == j
            continue
        else
            valid_index(1,pp) = j;
            select_matrix(1,pp) = train_target(i,:) * G * train_target(j,:)'; 
            pp = pp + 1;
        end
    end
    select_matrix = mapminmax(select_matrix,0,1);
    index_positive = find(select_matrix >= threshold);
    index_negative = find(select_matrix < threshold);
    candi_positive = select_matrix(1,index_positive);
    candi_negative = select_matrix(1,index_negative);
    [~,sort_positive] = sort(candi_positive,2);
    [~,sort_negative] = sort(candi_negative,2);
    num_positive_index = size(index_positive,2);
    num_negative_index = size(index_negative,2);

    if num_positive_index < num_similar
        diff_num = num_similar - num_positive_index;
        index_positive = [index_negative(sort_negative(1,end - diff_num + 1:end)),index_positive];
    end
    if num_negative_index < num_dissimilar
        diff_num = num_similar - num_negative_index;
        index_negative = [index_negative, index_positive(sort_positive(1,1:diff_num))];
    end

    valid_index_postive  = valid_index(1,index_positive);
    valid_index_negative = valid_index(1,index_negative);
    S_idex = knnsearch(train_data(valid_index_postive,:),train_data(i,:),'K',num_similar);
    R_idex = knnsearch(train_data(valid_index_negative,:),train_data(i,:),'K',num_dissimilar);
    Sneighbor(i,:) = valid_index_postive(S_idex);
    Rneighbor(i,:) = valid_index_negative(R_idex);
end

I = ones(num_label,num_label);
Dlmargin = I - (alpha * Anorm + beta * Bcorr);

%% optimization method :FISTA
[w] = COMMU_FISTA(train_data,train_target,Sneighbor,Rneighbor,Dlmargin,B,w,C);
M = eye(num_dim,num_dim);
for i=1:num_label
    M = M + w(i)*B(:,i)*B(:,i)';
end
end