function [prediction,beta,b]=MLUE(train_data,test_data,lambda1,lambda2,IFB,beta0,b0)
% Solve chain strcture L21 regularized multi-task learning with uncertainty estimation via EM algorithm
% [prediction,beta,b]=MLUE(train_data,test_data,lambda1,lambda2,IFB,beta0,b0)
% input:
% train_data: a struct, a collection of petition information from labeled petition sets.
% train_data has the following components:
% train_data.data: n * m labeled petition data matrix. -100 means this value is missing.
% train_data.index: n * 1 vector. The mapping from labeled petition data to labeled petition sets.
% train_data.label: n * 1 vector. Labels of labeled peititons. 1 means this petition is
% victorious, -1 means this petition fails.
% train_data.taskno: n * 1 vector. the task number assigned for multi-task learning according to
% missing patterns.
% test_data: a struct, a collection of petition information from unlabeled petition sets.
% test_data has the following components:
% test_data.data: k * m petition data matrix. -100 means this value is missing.
% test_data.index: k * 1 vector. The mapping from unlabeled petition data to unlabeled petition sets.
% test_data.taskno: k * 1 vector. The task number assigned for multi-task learning according to
% missing patterns.
% lambda1 and lambda2: two tuning parameters.
% IFB: 1*L cell, the partitions of increasing feature blocks.
% beta0 (optional): m*L matrix, initialization of beta.
% b0 (optional): 1*L vector, initialization of b.
% output:
% prediction: k *1 vector, the predicted values by the MLUE model.
% beta : m*L matrix, the cofficients for parameter.
% b: 1*L vector, the intercepts for L tasks.
% where:
% n: number of labeled petitions observations.
% m: number of features.
% k: number of unlabeled petitions observations.
% L: number of tasks.
max_iter=20;
fea_num=size(train_data.data,2);
max_taskno=max(train_data.taskno);
% if beta and b are included in the input, then use them as
% initializations. Otherwises, initialize them with random variables.
if exist('beta0','var')
    beta = beta0;
else
    rng('default');
    beta = 2*rand(fea_num,max_taskno)-1;
    for i=1:max_taskno
        unav_fea=setdiff(1:fea_num,IFB{i});
        beta(unav_fea,i)=0;
    end
end
if exist('b0','var')
    b = b0;
else
    rng('default');
    b = 2*rand(1,max_taskno)-1;    
end
Y=zeros(size(test_data.data,1),1);
% EM algorithm
for i=1:max_iter
    Y_last=Y;
    % update Y
    Y=update_Y(train_data,test_data,beta,b,lambda2,IFB);
    % update beta
    [beta,b]=update_beta(train_data,test_data,Y,beta,b,lambda1,IFB);
    % termination check
    if(all(Y_last==Y))
        break;
    end
end
% make predictions
prediction=zeros(size(test_data.data,1),1);
for i=1:max_taskno
prediction(test_data.taskno==i)=test_data.data(test_data.taskno==i,IFB{i})*beta(IFB{i},i)+b(i);
end

end