function [sep_measure,peak_order,transformation,decision_threshold] = cw_fisher(data,labels,num_genes)

dimPC = 3;
number = size (data, 1);
dimgene = size(data, 2);
dimsel = dimgene;
dimPeak = num_genes;

[id1,temp] = find(labels(:,1) == 1);
[id2,temp] = find(labels(:,1) == 2);
data1 = data(id1, :);
data2 = data(id2, :);




index=1;
vals=[];
for i=1:size(data1,1),
      if(length(vals) == 0 | length(find(vals(:,1) == labels(i,2))) == 0) 
        vals = [vals; labels(i,2),1];
        index=index+1;
      else
          [a,b]=find(vals(:,1) == labels(i,2));
          vals(a,2) = vals(a,2)+1;
      end
end
num_clusters(1)=index-1;
cluster1_num=vals(:,2)';


index=1;
vals=[];
offset=size(data1,1);
for i=1:size(data2,1),
      if(length(vals) == 0 |  length(find(vals(:,2) == labels(i+offset,2))) == 0) 
        vals = [vals; index, labels(i+offset,2),1];
        index=index+1;
      else
          [a,b]=find(vals(:,2) == labels(i+offset,2));
          vals(a,3) = vals(a,3)+1;

      end
end
num_clusters(2)=index-1;
cluster2_num=vals(:,3)';

pi_1 = cluster1_num/sum(cluster1_num);
pi_2 = cluster2_num/sum(cluster2_num);

num_samples(1)=size(data1,1);
num_samples(2)=size(data1,1);

%prescreening
M_inv_total = zeros(1,size(data1,2));
for k=1:size(data,2),
    [h1, mu1, sigmas1, P1] = em(num_clusters(1), data(1:num_samples(1),k));
    [h2, mu2, sigmas2, P2] = em(num_clusters(2), data(num_samples(1)+1:size(data,1),k));
    
    [new_labels,num_clusters,cluster1_num,cluster2_num] = get_cluster_labels(h1,h2);
    new_labels = [labels(:,1),new_labels];
    new_pi_1 = cluster1_num/sum(cluster1_num);
    new_pi_2 = cluster2_num/sum(cluster2_num);

  for i=1:num_clusters(1),
    [ind1,ind2] = find(new_labels(1:num_samples(1),2)== i);
    data_i=data1(ind1,:);
    Sw_i = cov(data_i);
    for j=1:num_clusters(2),
        [ind1,ind2] = find(new_labels(num_samples(1)+1:size(new_labels,1),2)== j);
        data_j=data2(ind1,:);

        stdcluster = (size(data_i(:,k), 1) - 1) * std(data_i(:,k), 0, 1).^2 + (size(data_j(:,k), 1) - 1) * std(data_j(:,k), 0, 1).^2;
        distance = (mean(data_i(:,k), 1) - mean(data_j(:,k), 1))^2;

        M_ij(k) = distance./stdcluster;
        
        M_inv_total(k) = M_inv_total(k) + (1./M_ij(k));
%        M_inv_total(k) = M_inv_total(k) + new_pi_1(i)*new_pi_2(j)*(1./M_ij(k));
    end 
end
end

 [temp, Peakorder] = sort(M_inv_total);
 Peakdata1 = data1(:, Peakorder);
 Peakdata2 = data2(:, Peakorder);
 Peakorder
 data1=Peakdata1(:,1:num_genes);
 data2=Peakdata2(:,1:num_genes);
 

%fisher discriminant calculation
M_total = 0;
for i=1:num_clusters(1),
    [ind1,ind2] = find(labels(1:num_samples(1),2)== i);
    data_i=data1(ind1,:);
    Sw_i = cov(data_i);
   for j=1:num_clusters(2),
        [ind1,ind2] = find(labels(num_samples(1)+1:size(labels,1),2)== j);
        data_j=data2(ind1,:);

        Sw_j = cov(data_j);
        
        Sw = (size(data_i,1)*Sw_i + size(data_j,1)*Sw_j)/(size(data_i,1) + size(data_j,1));
        M_ij = inv(Sw)*(mean(data_i, 1) - mean(data_j, 1))' * (mean(data_i, 1) - mean(data_j, 1));
        
        if(length(M_total)==1) M_total=zeros(size(M_ij,1),size(M_ij,2));
        end
         delta_ij = sqrt((mean(data_i, 1) - mean(data_j, 1)) * inv(Sw) * (mean(data_i, 1) - mean(data_j, 1))');
%          delta_ij = sqrt((mean(data_i, 1) - mean(data_j, 1)) * (mean(data_i, 1) - mean(data_j, 1))');
          M_total = M_total + M_ij*(erf(delta_ij/(2*sqrt(2)))/(2*delta_ij^2));
      end
end
    
[V, D] = eig(M_total);
[temp, PCorder] = sort(abs(diag(D)));
PCdata1 = data1 * V(:, PCorder(dimPeak - dimPC + 1 : dimPeak));
PCdata2 = data2 * V(:, PCorder(dimPeak - dimPC + 1 : dimPeak));

if(~isreal(V)) V;
end

IOdisplay(number, PCdata1, PCdata2, labels(:,1));


class1_mean = sum(PCdata1(:,3))/size(PCdata1,1);
class2_mean = sum(PCdata2(:,3))/size(PCdata2,1);
class1_std = std(PCdata1(:,3));
class2_std = std(PCdata2(:,3));
decision_threshold.val = (class1_mean*class2_std+class2_mean*class1_std)/(class1_std+class2_std);
%decision_threshold.val = (class1_mean+class2_mean)*.5;
if(class1_mean > decision_threshold.val) decision_threshold.label =1;
else decision_threshold.label=2;
end

% sep_measure = (mean(data1, 1) - mean(data2, 1)) * inv(Sw) * (mean(Peakdata1, 1) - mean(Peakdata2, 1))';
sep_measure=0;
% peak_order = Peakorder(dimsel - dimPeak + 1 : dimsel);
peak_order = [1:dimPeak];
transformation = V(:, PCorder(dimPeak - dimPC + 1 : dimPeak));