Data = importdata('gse6764.txt');
[nrow,ncol] = size(Data.data);
nordata = zeros(nrow,ncol); % normalized data
Data.data = log2(Data.data);
for i = 1:nrow
    for j = 1:ncol
        nordata(i,j) = (Data.data(i,j)-mean(Data.data(i,:)))/std(Data.data(i,:),1);
    end
end

file = textread('subnetwork6764.txt','%s','delimiter','\n');  %%% subnetwork6764.txt contains all the features (modules)
mdata = zeros(length(file),ncol);
for i = 1:length(file)
    words = strread(file{i},'%s','delimiter','\t');
    for j = 1:ncol
        for k = 1:length(words)
            if any(strcmp(words(k),Data.textdata(1:nrow+1)))==1
                ind = find(strcmp(words(k),Data.textdata(1:nrow+1))==1);
                mdata(i,j) = mdata(i,j)+nordata(ind-1,j);
            end
        end
    end
    mdata(i,j) = mdata(i,j)/sqrt(length(words));
end

st1 = 10; %% the number of the first-class samples
st2 = 13;
st3 = 17;
st4 = 18;
st5 = 17;

%%%%%%%%%%%%%%%%%%%%%%%%%% the center matrix of five-class samples
alpha=zeros(length(file),5);   
for i=1:length(file)
    alpha(i,1)=mean(mdata(i,1:st1));
    alpha(i,2)=mean(mdata(i,st1+1:st1+st2));
    alpha(i,3)=mean(mdata(i,st1+st2+1:st1+st2+st3));
    alpha(i,4)=mean(mdata(i,st1+st2+st3+1:st1+st2+st3+st4));
    alpha(i,5)=mean(mdata(i,st1+st2+st3+st4+1:st1+st2+st3+st4+st5));
end

diff = zeros(4*st1+4*st2+4*st3+4*st4+4*st5+2,length(file));
for j=1:length(file)
    for i=1:st1 % center1
        diff(i,j) = (mdata(j,i)-alpha(j,1))^2-(mdata(j,i)-alpha(j,2))^2;
        diff(i+st1,j) = (mdata(j,i)-alpha(j,1))^2-(mdata(j,i)-alpha(j,3))^2;
        diff(i+2*st1,j) = (mdata(j,i)-alpha(j,1))^2-(mdata(j,i)-alpha(j,4))^2;
        diff(i+3*st1,j) = (mdata(j,i)-alpha(j,1))^2-(mdata(j,i)-alpha(j,5))^2;
    end
    for i=1:st2 % center2
        diff(i+4*st1,j) = (mdata(j,i+st1)-alpha(j,2))^2-(mdata(j,i+st1)-alpha(j,1))^2;
        diff(i+4*st1+st2,j) = (mdata(j,i+st1)-alpha(j,2))^2-(mdata(j,i+st1)-alpha(j,3))^2;
        diff(i+4*st1+2*st2,j) = (mdata(j,i+st1)-alpha(j,2))^2-(mdata(j,i+st1)-alpha(j,4))^2;
        diff(i+4*st1+3*st2,j) = (mdata(j,i+st1)-alpha(j,2))^2-(mdata(j,i+st1)-alpha(j,5))^2;
    end
    for i=1:st3 % center3
        diff(i+4*st1+4*st2,j) = (mdata(j,i+st1+st2)-alpha(j,3))^2-(mdata(j,i+st1+st2)-alpha(j,1))^2;
        diff(i+4*st1+4*st2+st3,j) = (mdata(j,i+st1+st2)-alpha(j,3))^2-(mdata(j,i+st1+st2)-alpha(j,2))^2;
        diff(i+4*st1+4*st2+2*st3,j) = (mdata(j,i+st1+st2)-alpha(j,3))^2-(mdata(j,i+st1+st2)-alpha(j,4))^2;
        diff(i+4*st1+4*st2+3*st3,j) = (mdata(j,i+st1+st2)-alpha(j,3))^2-(mdata(j,i+st1+st2)-alpha(j,5))^2;
    end
    for i=1:st4 % center4
        diff(i+4*st1+4*st2+4*st3,j) = (mdata(j,i+st1+st2+st3)-alpha(j,4))^2-(mdata(j,i+st1+st2+st3)-alpha(j,1))^2;
        diff(i+4*st1+4*st2+4*st3+st4,j) = (mdata(j,i+st1+st2+st3)-alpha(j,4))^2-(mdata(j,i+st1+st2+st3)-alpha(j,2))^2;
        diff(i+4*st1+4*st2+4*st3+2*st4,j) = (mdata(j,i+st1+st2+st3)-alpha(j,4))^2-(mdata(j,i+st1+st2+st3)-alpha(j,3))^2;
        diff(i+4*st1+4*st2+4*st3+3*st4,j) = (mdata(j,i+st1+st2+st3)-alpha(j,4))^2-(mdata(j,i+st1+st2+st3)-alpha(j,5))^2;
    end
    for i=1:st5 % center5
        diff(i+4*st1+4*st2+4*st3+4*st4,j) = (mdata(j,i+st1+st2+st3+st4)-alpha(j,5))^2-(mdata(j,i+st1+st2+st3+st4)-alpha(j,1))^2;
        diff(i+4*st1+4*st2+4*st3+4*st4+st5,j) = (mdata(j,i+st1+st2+st3+st4)-alpha(j,5))^2-(mdata(j,i+st1+st2+st3+st4)-alpha(j,2))^2;
        diff(i+4*st1+4*st2+4*st3+4*st4+2*st5,j) = (mdata(j,i+st1+st2+st3+st4)-alpha(j,5))^2-(mdata(j,i+st1+st2+st3+st4)-alpha(j,3))^2;
        diff(i+4*st1+4*st2+4*st3+4*st4+3*st5,j) = (mdata(j,i+st1+st2+st3+st4)-alpha(j,5))^2-(mdata(j,i+st1+st2+st3+st4)-alpha(j,4))^2;
    end
end
diff(4*st1+4*st2+4*st3+4*st4+4*st5+1,:) = -ones(1,length(file));
diff(4*st1+4*st2+4*st3+4*st4+4*st5+2,:) = ones(1,length(file));
b = zeros(4*st1+4*st2+4*st3+4*st4+4*st5+2,1);
b(4*st1+4*st2+4*st3+4*st4+4*st5+1,1) = -1;
lb = zeros(length(file),1);
ub = ones(length(file),1);

lamda = 0.01:0.1:1; % lamda should be given.
lowerno = 8;
upno = 30;
pow = zeros(length(lamda),upno-lowerno+1);
for no = lowerno:upno
    b(4*st1+4*st2+4*st3+4*st4+4*st5+2,1) = no;
    for i =1:length(lamda)
        f = zeros(length(file),1); % sumj(1+lamda*sum(diff(:,j)))xj
        for j = 1:length(file)
            f(j,1) = 1+lamda(i)*(sum(diff(1:4*st1+4*st2+4*st3+4*st4+4*st5,j)));
        end
        [x,fval] = linprog(f,diff,b,[],[],lb,ub);
        for k=1:length(file)
            if x(k)>0.9
                x(k)=1;
            else
                x(k)=0;
            end
        end
        idiff = diff(1:4*st1+4*st2+4*st3+4*st4+4*st5,:)*x;
        num = 0;
        for ll=1:st1 % center1
            if idiff(ll)>0 || idiff(ll+st1)>0 || idiff(ll+2*st1)>0 || idiff(ll+3*st1)>0
                num = num+1;
            end
        end
        for ll=1:st2 % center1
            if idiff(ll+4*st1)>0 || idiff(ll+4*st1+st2)>0 || idiff(ll+4*st1+2*st2)>0 || idiff(ll+4*st1+3*st3)>0
                num = num+1;
            end
        end
        for ll=1:st3 % center1
            if idiff(ll+4*st1+4*st2)>0 || idiff(ll+4*st1+4*st2+st3)>0 || idiff(ll+4*st1+4*st2+2*st3)>0 || idiff(ll+4*st1+4*st2+3*st3)>0
                num = num+1;
            end
        end
        for ll=1:st4 % center1
            if idiff(ll+4*st1+4*st2+4*st3)>0 || idiff(ll+4*st1+4*st2+4*st3+st4)>0 || idiff(ll+4*st1+4*st2+4*st3+2*st4)>0 || idiff(ll+4*st1+4*st2+4*st3+3*st4)>0 
                num = num+1;
            end
        end
        for ll=1:st5 % center1
            if idiff(ll+4*st1+4*st2+4*st3+4*st4)>0 || idiff(ll+4*st1+4*st2+4*st3+4*st4+st5)>0 || idiff(ll+4*st1+4*st2+4*st3+4*st4+2*st5)>0 || idiff(ll+4*st1+4*st2+4*st3+4*st4+3*st5)>0 
                num = num+1;
            end
        end
        pow(i,no-lowerno+1) = num;
    end
end
indno = find(min(pow(1,:))==pow(1,:));
b(4*st1+4*st2+4*st3+4*st4+4*st5+2,1) = indno(1)+8-1;
f = zeros(length(file),1); % sumj(1+lamda*sum(diff(:,j)))xj
for j = 1:length(file)
    f(j,1) = 1+lamda(i)*(sum(diff(1:4*st1+4*st2+4*st3+4*st4+4*st5,j)));
end
[x,fval] = linprog(f,diff,b,[],[],lb,ub);
for k=1:length(file)
    if x(k)>0.9
        x(k)=1; % it means that the k-th module is selected.
    else
        x(k)=0; % it means that the k-th module is not selected.
    end
end

indmod = find(x==1);
score = zeros(length(indmod),ncol);
for i = 1:length(indmod)
    for j = 1:ncol
        score(i,j) = mdata(indmod(i),j);
    end
end
fidscore = fopen('HCCidemodscore.txt','w');
for i = 1:length(indmod)
    for j = 1:ncol
        fprintf(fidscore,'%f\t',score(i,j));
    end
    fprintf(fidscore,'\n');
end

