function [test_auc,cv_auc] = pooling(C, alpha, maxiter, data, labels, testdata, testlabels, cvdata, cvlabels,ktype,kdeg,par)
gamma1=1;
num_tasks = size(data,2);
ctot = 0;
ttot = 0;
tot = 0;
for i=1:num_tasks
  ctot = ctot + size(cvdata{i},1);
  ttot = ttot + size(testdata{i},1);
  tot = tot + size(data{i},1);
end

lambdas = zeros(tot,1);

Y = [];
T = [];

CVY = []; %CV Y
CVT = []; %CV test

TEY = []; %CV Y
TET = []; %CV test

param.MSK_DPAR_INTPNT_CO_TOL_PFEAS = 1.0e-15;

% Dual feasibility tolerance for the dual solution
param.MSK_DPAR_INTPNT_CO_TOL_DFEAS = 1.0e-15;

% Relative primal-dual gap tolerance.
param.MSK_DPAR_INTPNT_CO_TOL_REL_GAP = 1.0e-15;
param.MSK_IPAR_INTPNT_NUM_THREADS = 2;


for m = 1:num_tasks % for each pool
    Y = [Y; labels{m}];
    T = [T; data{m}];

    CVY = [ CVY; cvlabels{m} ];
    CVT = [ CVT; cvdata{m} ];

    TEY = [ TEY; testlabels{m} ]; 
    TET = [ TET; testdata{m} ]; 
end

X = par*makekernel(T,T,ktype,kdeg).*(Y*Y');
CX= par*makekernel(T,CVT,ktype,kdeg).*(Y*ones(1,length(CVY)));
TX =par*makekernel(T,TET,ktype,kdeg).*(Y*ones(1,length(TEY)));


%%% compute H_m, CH_m and TH_m

curr = 0; % offset
task = 1;
for m = 1:1:size(data,2) % for each pool
    H{m} = zeros(tot); % empty tot x tot matrix
    for t = 1:size(data{m}, 1)
        for tau = 1:size(data{m}, 1)
            H{m}(curr+t, curr+tau) = (1-par)*num_tasks*labels{m}(t)*labels{m}(tau)*makekernel(data{m}(t,:),data{m}(tau,:),ktype,kdeg); %y_t * y_tau * x_t * x_tau
        end
    end
    curr = curr+size(data{m},1); % increase offset
    CH{m} =(1-par)*num_tasks*makekernel(data{m},cvdata{m},ktype,kdeg).*(labels{m}*ones(1,length(cvlabels{m})));
    TH{m} =(1-par)*num_tasks*makekernel(data{m},testdata{m},ktype,kdeg).*(labels{m}*ones(1,length(testlabels{m}))) ;
end

going = 1;
it = 0;

% iterate solution
fval = inf;
while(going && it < maxiter)
    % compute s and g

    for m = 1:1:size(data,2)
        if (alpha == inf)
            temp = inf;
        else
            temp = alpha*exp(-(lambdas'*H{m}*lambdas)/2);
        end    
        s(m) = 1/(temp+1);
        if (temp==1) 
            g(m)=0;
        else
            g(m) = tanh(0.5*log(temp))/(2*log(temp));
        end    
    end

    % linear term

    
    prob.c = -gamma1*ones(1,tot);
    for m = 1:1:size(data,2)
        prob.c = prob.c + s(m)*lambdas'*H{m} - lambdas'*H{m} - g(m)*lambdas'*H{m}*lambdas*lambdas'*H{m};
    end

    
    % quadratic term

    
    Q = X;
    for m = 1:1:size(data,2) % for each pool
        Q = Q + g(m) * H{m} * lambdas * lambdas' * H{m} + H{m};
    end
    
    Q=(Q+Q')/2+1e-6*eye(tot); %avoid numerical problems

    prob.blx = zeros(tot,1);
    prob.bux = C*ones(tot,1);


    %build A
    A = zeros(size(data,2), tot);
    offset = 0;
    for m = 1:size(data,2) % for each class
        for t = 1:size(data{m},1) % for each point in m
            A(m,t+offset) = labels{m}(t);
        end
        offset = offset + size(data{m},1);
    end

    prob.a = sparse(A);     
    prob.buc = zeros(size(data,2),1);
    prob.blc = zeros(size(data,2),1);
    old = lambdas;

    
    % maximize: -H and -f
    oldfval = fval;
    [prob.qosubi,prob.qosubj,prob.qoval] = find(tril(sparse(double(Q))));
    [r,res]=mosekopt('minimize echo(0)',prob,param);
    lambdas = res.sol.itr.xx(1:tot);
    bias =  res.sol.itr.suc - res.sol.itr.slc;
    fval = 0.5*lambdas'*Q*lambdas + prob.c*lambdas;
    %assert( fval < oldfval ); 
    %fprintf('%f %f\n',fval,norm(old-lambdas)/norm(lambdas));

    if (norm(old-lambdas)/norm(lambdas) < 1e-6)
        going = 0;
    end
    it=it+1;
end

for m = 1:1:size(data,2)
    temp = alpha*exp(-(lambdas'*H{m}*lambdas)/2);
    s(m) = 1/(temp+1);
    fprintf('%f ',s(m));
    %fprintf('%f   %f      %f\n',alpha,exp(-(lambdas'*H{m}*lambdas)/2),s(m));
end
fprintf('\n');


%%Calculate the test and cv predictions (without the bias term)
test_pred = [];
cv_pred   = [];
offset = 0;
coffset = 0;
toffset = 0;
test_labels_all = [];
cv_labels_all = [];
for m=1:size(data,2)
  test_labels_all = [ test_labels_all; testlabels{m} ];
  cv_labels_all = [ cv_labels_all; cvlabels{m} ];

  ind = offset+1:offset+size(data{m},1);
  cvind = coffset+1:coffset+size(cvdata{m},1);
  testind = toffset+1:toffset+size(testdata{m},1);
 
  temp = s(m)*(CH{m})'*lambdas(ind);
  temp = temp + CX(:,cvind)'*lambdas ;
  cv_pred = [ cv_pred ; temp + bias(m) ];

  temp = s(m)*(TH{m})'*lambdas(ind);
  temp = temp + TX(:,testind)'*lambdas;
  test_pred = [ test_pred; temp + bias(m) ];

  offset = offset + size(data{m},1); 
  coffset = coffset + size(cvdata{m},1); 
  toffset = toffset + size(testdata{m},1); 
end


%% compute the AUC
test_auc = roc_auc(test_pred, test_labels_all);
cv_auc   = roc_auc(cv_pred,cv_labels_all);
