function demo_ap_cluster(N, num_pref)
    
    addpath(genpath('~/code/MrCat-dev'))
    
    mu1 = [7 8];
    sigma1 = [1 .1; .1 .3];
    
    mu2 = [3 3];
    sigma2 = [.4 0; 0 .5];
    
    mu3 = [8 3];
    sigma3 = [1 .1; .1 .3];
       
%     X = [mvnrnd(mu1, sigma1, round(N/3)); 
%            mvnrnd(mu2, sigma2, round((N - round(N/3))/2)); 
%            mvnrnd(mu3, sigma3, (N-(round((N - round(N/3))/2)+round(N/3))))];

%   X = generateData(1, 0.5, clusters, 15, 15, 5, 1, 2, N);

X = [randn(N,2)*0.2+ones(N,2)+2;
     randn(N,2)*0.3-ones(N,2)-2;
     randn(N,2)*0.4+ones(N,2)+4;
     randn(N,2)*0.2-ones(N,2)-4;
     randn(N,2)*0.3+ones(N,2)+6;
     randn(N,2)*0.4-ones(N,2)-6;
     randn(N,2)*0.2+ones(N,2)+8;
     randn(N,2)*0.3-ones(N,2)-8];

    % Visualise raw data
    figure;
    g=gramm('x', X(:,1), 'y',X(:,2));
    g.geom_point();
    g.set_title(['simulated data of ' num2str(N) ' data points']);
    g.draw();
    
    x = X;
   % x = rand(N,2);
    M = N*N-N;
    s = zeros(M,3);
    j = 1;

    for i = 1:N
        for k = [1:i-1,i+1:N]
            s(j,1) = i;
            s(j,2) = k;
            s(j,3) = -sum((x(i,:)-x(k,:)).^2);
            j=j+1;
        end;
    end;

    fprintf('Calculating range of preference values \n');
    [pmin,pmax] = preferenceRange(s,'exact');
    pref_values = linspace(pmin,pmax,num_pref);    

    num_clusters = [];
    cluster_arrangement = [];
    net_sim = []; dp_sim = []; ex_pref = [];

    for ipref = 1: length(pref_values)
        if mod(ipref, 10)==1
            fprintf('Running loop: %d (of %d)\n', ipref, num_pref)
            fprintf('Preference value: %d \n', pref_values(ipref))
        end

        %    [idx,netsim,dpsim,expref] = apcluster(data_CC_ap, pref_values(ipref), 'plot','nonoise');
        [idx,netsim,dpsim,expref] = apcluster(s,pref_values(ipref),'nonoise');

        net_sim = [net_sim; netsim];
        dp_sim = [dp_sim; dpsim];
        ex_pref = [ex_pref; expref];
        cluster_arrangement = [cluster_arrangement; idx];
        this_solution = length(unique(idx));
        num_clusters = [num_clusters; this_solution];

    end
    pref_values = pref_values';    
    
    %%    
    % Examine specific preference-value-cases 
    
    % median(s)
    p_median = median(s(:,3));
    [idx,netsim,~,~] = apcluster(s,p_median,'nonoise');
    num_median_clusters = numel(unique(idx));
    net_sim = [net_sim; netsim];
    
    % min(s)
    p_min = min(s(:,3));
    [idx,netsim,~,~] = apcluster(s,p_min,'nonoise');
    num_min_clusters = numel(unique(idx));
    net_sim = [net_sim; netsim];
    
    % min(s) - (max(s) - min(s))
    p_strict = min(s(:,3))-(max(s(:,3))-min(s(:,3)));
    [idx,netsim,~,~] = apcluster(s,p_strict,'nonoise');
    num_strict_clusters = numel(unique(idx));
    net_sim = [net_sim; netsim];
    
    %%
    data = [];

    % Group preferences for plotting in different colours 
    preference_type = ones(1,length(pref_values))';

    data.pref_values    = vertcat(pref_values, p_median, p_min, p_strict);
    data.num_clusters   = vertcat(num_clusters, num_median_clusters, num_min_clusters, num_strict_clusters);
    data.preference_type      = vertcat(preference_type, 2, 3, 4);
    data.net_sim        = net_sim;

    for i=1:length(data.preference_type)
        if data.preference_type(i) == 1
            data.lable{i} = 'a pref value';
        elseif data.preference_type(i) == 2
            data.lable{i} = 'median(S)';
        elseif data.preference_type(i) == 3
            data.lable{i} = 'min(S)';
        elseif data.preference_type(i) == 4
            data.lable{i} = 'min(S) - (max(S) - min(S))';            
        end
    end
    
    save all_pref_test
    save('data', 'data')
    
    figure;
    g=gramm('x', data.pref_values, 'y', data.num_clusters, 'color', data.lable);
    g.geom_jitter();
    g.set_names('x', 'preference values', 'y','number of clusters','color','what preference value?');
    y_limit = unique(data.num_clusters);
    g.axe_property('YLim',[0 y_limit(end-1)+2]);
    % g.axe_property('XLim',[data.pref_values(1) data.pref_values(end-1)]);
    g.set_title('preference value controls number of clusters');
    g.draw();
    
    figure;
    g=gramm('x', data.pref_values, 'y', data.net_sim, 'color', data.lable);
    g.geom_point();
    g.set_names('x', 'preference values', 'y','net similarity','color','what preference value?');    
    g.set_title('net similatiry - optimization parameter for clustering');
    g.draw();

end
