function [data,labels] = generate_gmm_data(parameters)

dim = parameters.dim;
num_samples(1) = parameters.num_samples(1);
num_samples(2) = parameters.num_samples(2);
num_clusters(1) = parameters.num_clusters(1);
num_clusters(2) = parameters.num_clusters(2);
num_indicators = parameters.num_indicators;
indicator_strength = parameters.indicator_strength;
indicator_var = parameters.indicator_var;
data_range = parameters.range;
max_var = parameters.max_var;

num1_left = num_samples(1);
cluster1_num=zeros(1,num_clusters(1));
while(num1_left > 0)
    for i=1:num_clusters(1),
        temp = ceil(rand*(num1_left/num_clusters(1)));
        if(temp > num1_left) temp = num1_left;
        end
        cluster1_num(i) = cluster1_num(i) + temp;
        num1_left = num1_left - temp;
        if(num1_left == 0) i = num_clusters(1); break;
        end
    end
end

num2_left = num_samples(2);
cluster2_num=zeros(1,num_clusters(2));
while(num2_left > 0)
    for i=1:num_clusters(2),
        temp = ceil(rand*(num2_left/num_clusters(2)));
        if(temp > num2_left) temp = num2_left;
        end
        cluster2_num(i) = cluster2_num(i) + temp;
        num2_left = num2_left - temp;
        if(num2_left == 0) i = num_clusters(2); break;
        end
    end
end

scale = .25;
class1_centroid = data_range(1) + randn(1,num_indicators)*(data_range(2)-data_range(1));
class2_centroid = data_range(1) + randn(1,num_indicators)*(data_range(2)-data_range(1));
dist = norm(class1_centroid-class2_centroid);
while(dist < indicator_strength | dist > 2*indicator_strength)
    if(dist < indicator_strength) class2_centroid = class2_centroid + (class2_centroid - class1_centroid)*scale;
    else class2_centroid = class2_centroid - (class2_centroid - class1_centroid)*scale;
    end
    dist = norm(class1_centroid-class2_centroid);
end


    class1=[];        
for i=1:num_clusters(1),
    mean1(i,:) = mvnrnd(class1_centroid,diag(zeros(1,num_indicators) + max_var*5),1);
    class1 = [class1; mvnrnd(mean1(i,:),diag(zeros(1,num_indicators) + rand*max_var),cluster1_num(i))];
end

class2=[];
for i=1:num_clusters(2),
    mean_seps = normrnd(indicator_strength,indicator_var,1,num_clusters(1));
    dist = zeros(1,num_clusters(1));
%     while(sum(dist < .5*mean_seps | dist > 4*mean_seps) > 0)
        mean2(i,:) = mvnrnd(class2_centroid,diag(zeros(1,num_indicators) + max_var*5),1);
%         for j=1:num_clusters(1),
%             dist(j) = norm(mean2(i,:) - mean1(j,:));
%         end
%     end
    class2 = [class2; mvnrnd(mean2(i,:),diag(zeros(1,num_indicators) + rand*max_var),cluster2_num(i))];
end

            
            
%     mean1(i) = data_range(1) + randn(1,num_indicators)*(data_range(2)-data_range(1));
%     class1 = [class1, mvnrnd(mean1,diag(zeros(1,num_indicators) + rand*(data_range(2)-data_range(1))/5),cluster1_num(i))];
% end



% for i=1:num_indicators,
%     mean_sep = normrnd(indicator_strength,indicator_var);
%     mean1 = data_range(1) + rand*(data_range(2)-data_range(1));
%     mean2 = mean1 + mean_sep;
%     class1(:,i) = normrnd(mean1,rand*(data_range(2)-data_range(1))/5,num_samples(1),1);
%     class2(:,i) = normrnd(mean2,rand*(data_range(2)-data_range(1))/5,num_samples(2),1);
% end


for i=1:(dim-num_indicators),
    rand_mean = data_range(1) + rand*(data_range(2)-data_range(1));
    rand_var = rand*(data_range(2)-data_range(1))/5;
    temp = normrnd(rand_mean,rand_var,num_samples(1)+num_samples(2),1);
    class1(:,i+num_indicators) = temp(1:num_samples(1));
    class2(:,i+num_indicators) = temp(num_samples(1)+1:length(temp));
end


% temp=randn(1,dim);
% [v,i]=sort(temp);
i=[1:size(class1,2)];
data(1:num_samples(1),:) = class1(:,i);
data(num_samples(1)+1:sum(num_samples),:) = class2(:,i);


index=1;
for i=1:num_clusters(1),
   for j=1:cluster1_num(i),
       labels(index,:)=[1,i];
       index=index+1;
   end
end

for i=1:num_clusters(2),
   for j=1:cluster2_num(i),
       labels(index,:)=[2,i];
       index=index+1;
   end
end
