num_samples=[100,100];
num_clusters=[3,4];
num_genes=200;
num_indicators=100;
ind_strength = (10)*num_genes;

[data,labels] = gmm(num_samples,num_clusters,num_genes,num_indicators,ind_strength);

n=50;
lambda_thresh = .1;

[inf_genes,vals] = get_max_fisher_genes(data,labels(:,1),n);
lambda_avg = mean(vals);

reduced_data = data(:,inf_genes);
for k=1:size(reduced_data,2),
  for i=1:num_clusters(1),
    [ind1,ind2] = find(labels(1:num_samples(1),2)== i);
    data_i=data(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=data(num_samples(1)+ind1,:);

        stdcluster = (size(data_i, 1) - 1) * std(data_i(:,k), 0, 1).^2 + (size(data_j, 1) - 1) * std(data_j(:,k), 0, 1).^2;
        distance = (mean(data_i(:,k), 1) - mean(data_j(:,k), 1)).^2;
        lambda = distance/stdcluster;
        
        if(i==1 & j==1) min_lambda(k,:) = [i,j,lambda];
        elseif(lambda < min_lambda(k,3)) min_lambda(k,:) = [i,j,lambda];
        else ;    
        end
    end
  end
end

[ind1,ind2] = find(min_lambda(:,3) < lambda_thresh*lambda_avg);
