function [W pi P obj] = LogisticL2RegressionR1groupR21_different(X,Y,Y1,group_idx,r1,r2)
% X: each column is a data point
% Y,Y1: matrices for the two regressions. n by c1(c2) matrices.
%    for logistic regression, each row of Y denotes a class indicator vector.
%    if xi belongs to class k, then Y(i,k)=1 and Y(i,k1)=0(k1!=k)
% group_idx: group index of features
% r1, r2: regularization parameter

% W: output regression matrix for logistic regression
% P: output regression matrix for least squares regression

% solve the following problem:
% min_{V=[W P]} Logit Loss(W,X,Y) + ||X'P-Y1||^2 + r1*||V||_1group + r2*||V||_21


group_set = unique(group_idx);
group_num = length(group_set);

[dim,n] = size(X);
class_num = size(Y,2);
y = sum(Y,2)';
YY = ones(dim,1)*y;

class_num1 = size(Y1,2);


%W = 0.1*rand(dim,class_num); P = 0.1*rand(dim,class_num1);
%W = zeros(dim,class_num); P = zeros(dim,class_num1);
W = ones(dim,class_num); P = ones(dim,class_num1);
V = [W P];


dimall = dim*class_num;

XY = X*Y;
XY1 = X*Y1;
XX = X*X';

t = cputime;
%for iter = 1:150
for iter = 1:10

    for c = 1:group_num
        idx = find(group_idx==group_set(c));
        Vc = V(idx,:);
        di = sqrt(sum(Vc.*Vc,1)+eps);
        Vig(idx,:) = ones(length(idx),1)*di;
        ob(c) = sum(di);
    end;
    D1 = 0.5./Vig;
    Vi = sqrt(sum(V.*V,2)+eps);
    D2 = diag(0.5./(Vi));
    
    for i = 1:class_num1
        P(:,i) = (XX+r1*diag(D1(:,class_num+i))+r2*D2)\XY1(:,i);
    end;
    
    
    WX = W'*X;
    t1 = exp(WX);
    t2 = sum(t1);
    tt = ones(class_num,1)*t2;
    pi = t1./tt;
    
    te1 = WX'.*Y;
    te1 = sum(te1(:));
    te2 = log(tt').*Y;
    te2 = sum(te2(:));
    XPY = X'*P - Y1;
    obj(iter) = te2 - te1 + trace(XPY'*XPY) + r1*sum(ob) + r2*sum(Vi);
      
    d = zeros(dim,class_num);
    for k = 1:class_num
        temp = ones(dim,1)*pi(k,:);
        G = YY.*X.*temp;
        d(:,k) = sum(G,2);
    end;
    d = d - XY;
    Wd1 = D1(:,1:class_num).*W;
    Wd2 = D2*W;
    d = d(:) + 2*r1*Wd1(:) + 2*r2*Wd2(:);
    
    H = zeros(dimall);
    for k = 1:class_num
        kdix = (k-1)*dim+1:k*dim;
        for m = 1:class_num
            mdix = ((m-1)*dim+1):(m*dim);
            D = sparse(diag(y.*pi(k,:).*pi(m,:)));
            if k == m
                D11 = sparse(diag(y.*pi(k,:)));
                H(kdix, mdix) = X*D11*X' - X*D*X' + 2*r1*diag(D1(:,k)) + 2*r2*D2;
            else
                H(kdix, mdix) = -X*D*X';
            end;
        end;
    end;
    
    b = H\d;
    W = W(:) - b;
    W = reshape(W,dim,class_num);
    V = [W P]; 
	
	disp(['The ', num2str(iter, '%02d'), 'th iteration has been done by ', num2str(cputime - t, '%06.2f'), ' seconds, obj_val is ', num2str(obj(iter))]);
end;




