function [w_f, alpha, alpha_update, obj_hist,Y_infer,hamm_loss] = BCD_prox_gradient_CGM_real(prob,  model, param, max_iter)
% Master problem optimization w.r.t. the sequared hinge loss in the primal form
%
%Input:
%   prob:      containing some basics of the problem
%   model:
%   param:
% Written by Mingkui Tan. If any comments, please email tanmingkui@gmail.com
%%
if nargin < 4
    max_iter = 50;
end

%%
prob.norm_sq = sum(prob.X.^2,2);
lambda = param.lambda;
out_put_loss = param.out_put_loss;
update_close_form = param.update_close_form;
G = model.G; %
alpha = model.alpha;
num_cliques = model.num_cliques; %cliques of order higher than 1 

X = prob.X;
y = prob.y;
D = 0;
N_inst = prob.num_inst;
Dim = prob.Dim;
num_label = prob.num_label;
lambda_N_inst = lambda*N_inst;


total_num_cliques = num_label + num_cliques;
total_num_para = 2*num_label+1*num_cliques;

if num_cliques==0
    w_f = zeros(num_label*2, Dim);
    w_inst = w_f;
else
    w_f = model.w;
    w_inst = w_f;   
end
weight = ones(size(w_f,1),1);
weight(2*num_label+1:end,1) = 1;
%%
% SG = zeros(n_feat,1);
% normsq = zeros(1,num_cliques);
norm_sq  = max(prob.norm_sq);


f_obj = 0;
obj_hist = zeros(max_iter,1);
Hamm_loss = zeros(N_inst,1);

i = 1;
x_inst = X(i,:);
y_inst = y(i,:);


[y_infer, loss_inst, delta_phi, hamm_y_yj]= infer_most_vialated_feature_mapping_v1(prob, G, w_f, x_inst, y_inst);
S = lambda_N_inst*loss_inst/norm(x_inst)^2/total_num_cliques;
if S>0.1
    scale = max(2, S/0.1);
else
    scale = 4;
end

for iter = 1:max_iter
    
    rand_seq = randperm(N_inst);
    f_obj = 0;
%   alpha = zeros(N_inst,1);
    for i_inst = 1:N_inst
        
        % choose a point
        i = rand_seq(i_inst);
        x_inst = X(i,:);
        y_inst = y(i,:);

        % infer y_infer, and  compute_delta_phi_y defines specific feature functions
        [y_infer, loss_inst, delta_phi, hamm_y_yj]= infer_most_vialated_feature_mapping_v1(prob, G, w_f, x_inst, y_inst);
        Hamm_loss(i,1) = hamm_y_yj;
        Y_infer(i,:) = y_infer;
        f_obj = f_obj + loss_inst;
        
        %compute step size S, update alpha and w
        if update_close_form > 0
             S = lambda_N_inst*loss_inst/prob.norm_sq(i)/total_num_para/scale;
             alpha_temp = alpha(i);
             alpha(i) = min(1,max(alpha_temp + S,0));
             S = alpha(i) - alpha_temp;
             non_zero_delta_y = find(delta_phi~=0);
             w_f(non_zero_delta_y,:) = w_f(non_zero_delta_y,:) + S*kron(delta_phi(non_zero_delta_y,:),x_inst)/lambda_N_inst;     
        end
    end

        if out_put_loss == 1
            obj_hist(iter,1) =   lambda*sum(sum(w_f.*w_f))/2 + 1/N_inst*f_obj;
        end
        if iter>2
            if (obj_hist(iter-1,1) - obj_hist(iter,1))/obj_hist(iter-1,1) < 0.00001
                scale = 2*scale;
            end
        end
        if (iter>5)
            ratio = abs(obj_hist(iter,1)-obj_hist(iter-5,1))/abs(obj_hist(iter-5,1));
            if (ratio<0.001)
            break
            end
        end
    %do acceleration if necessary
end
obj_hist(iter:end) = [];
hamm_loss = 0;

%   update alpha_iy after optimization
for j = 1:N_inst
     x_inst = X(j,:);
     y_inst = y(j,:);
     [y_infer, loss_inst, delta_phi, hamm_y_yj]= infer_most_vialated_feature_mapping_v1(prob, G, w_f, x_inst, y_inst);
     hamm_loss = hamm_loss + hamm_y_yj;
     alpha(j) = min(loss_inst, num_label);
     Y_infer(j,:) = y_infer;
end
alpha_update = alpha;
hamm_loss = hamm_loss/N_inst/num_label;




