
X = XX;
Phi2 = Phi;
min_dst = sqrt(2.2)*1/ (sqrt(N)-1);
dd = pdist2(X, X);
dd = 10* tril(dd) + dd + 5*eye(N);
neigh = dd < min_dst;
ii = find(neigh > 0);
[i1, i2] = ind2sub(size(neigh),ii);
D = zeros(size(i1,1), 1);
Y2 =  Phi2 * alphas;


MSEY  = @(idx) mean(   (Eta(idx, :) - mean(Eta(idx, :)) ).^2, 'all');    

tot = zeros(niter- burnin, 2);
for iii=1:(niter- burnin)
    if mod(iii,10) == 0
        disp(iii);
    end
    [Evec2, Eval2] = eig(Omegas(:,:,iii), "vector");
    [Eval,ind] = sort(Eval2, 'descend');
    Evec = Evec2(:, ind);
    E1 = Evec * diag(sqrt(Eval));
    Eta = Phi2 * E1;
    D  = sum(abs(Eta(i1, :) - Eta(i2, :)), 2);  %vecnorm(Eta(i1, :) - Eta(i2, :), 2, 2);








    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  MST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%







    G = graph(i1, i2, D);
    [T,pred] = minspantree(G);
    %LWidths = 5*T.Edges.Weight/max(T.Edges.Weight);
    %plot(T, 'XData', XX(:,1), 'YData', XX(:,2), 'LineWidth', LWidths);











    min_clust_size = N*0.5/100;
    max_clust_size = N*2/100;
    [~, idx] = sort(T.Edges.Weight, 'descend');


    IDX = [];
    j_cur = 0;
    maxcnt = 2000;
    while size(IDX) < 150
        for j=(1+j_cur):(N-1)
            tmp_IDX = [IDX idx(j)];
            U  = T.rmedge(tmp_IDX);
            if maxcnt <= max_clust_size
                break;
            end
            mincnt  = min(   histc(U.conncomp,  unique(U.conncomp))    );
            %maxcnt  = max(   histc(U.conncomp,  unique(U.conncomp))    );
            if mincnt >= min_clust_size
                IDX = tmp_IDX;
                j_cur = j;
                %disp(j);
                break
            end
            if j >= N*90/100
                break;
            end
        end
        if j >= N*90/100
            break;
        end
        if maxcnt <= max_clust_size
                break;
        end
    end
    
    
    
    
    U = T.rmedge(IDX);
    uu = U.conncomp';
    MSEs = zeros(max(uu), 1); 
    for i=1:max(uu)
        idx1 = find(uu == i);
        MSEs(i) = MSEY(idx1);    
    end
    tot(iii, :)  = [size(IDX,2), sum(MSEs)];
end


[~, iii] = min(  tot(:,2) ./ tot(:,1) );
[Evec2, Eval2] = eig(Omegas(:,:,iii), "vector");
[Eval,ind] = sort(Eval2, 'descend');
Evec = Evec2(:, ind);
E1 = Evec * diag(sqrt(Eval));
Eta = Phi2 * E1;
D  = sum(abs(Eta(i1, :) - Eta(i2, :)), 2);  %vecnorm(Eta(i1, :) - Eta(i2, :), 2, 2);








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  MST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%







G = graph(i1, i2, D);
[T,pred] = minspantree(G);
%LWidths = 5*T.Edges.Weight/max(T.Edges.Weight);
%plot(T, 'XData', XX(:,1), 'YData', XX(:,2), 'LineWidth', LWidths);











min_clust_size = N*0.5/100;
max_clust_size = N*2/100;
[~, idx] = sort(T.Edges.Weight, 'descend');


IDX = [];
j_cur = 0;
maxcnt = 2000;
while size(IDX) < 150
    for j=(1+j_cur):(N-1)
        tmp_IDX = [IDX idx(j)];
        U  = T.rmedge(tmp_IDX);
        if maxcnt <= max_clust_size
            break;
        end
        mincnt  = min(   histc(U.conncomp,  unique(U.conncomp))    );
        %maxcnt  = max(   histc(U.conncomp,  unique(U.conncomp))    );
        if mincnt >= min_clust_size
            IDX = tmp_IDX;
            j_cur = j;
            %disp(j);
            break
        end
        if j >= N*90/100
            break;
        end
    end
    if j >= N*90/100
        break;
    end
    if maxcnt <= max_clust_size
            break;
    end
end




U = T.rmedge(IDX);
uu = U.conncomp';

Yhat = Y2(:,1);
aehat = zeros(N,1);

for j =1:size(unique(U.conncomp),2)
    Yhat(U.conncomp == j, 1) = mean(Y2(U.conncomp == j, iii), "all");
end

e  = abs(Yhat - YY);
ae = e;
for j =1:size(unique(U.conncomp),2)
    ae(U.conncomp == j, 1) = sqrt(mean(e(U.conncomp == j, 1).^2));
end








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



bottomY = min([YY, Yhat], [], 'all');
topY    = max([YY, Yhat], [], 'all');

bottome = min(e);
tope    = max(e);



%f = subplot(2,2);
%nexttile;
subplot(2,2,1);
scatter(X(:,1), X(:,2), 2, YY);
shading interp;
caxis manual
caxis([bottomY topY]);
title('Point Level Data');
colorbar;
%nexttile;
subplot(2,2,2);
scatter(X(:,1), X(:,2), 2, Yhat);
shading interp;
caxis manual
caxis([bottomY topY]);
title('Areal Level Data');
colorbar;
%nexttile;
subplot(2,2,3);
scatter(X(:,1), X(:,2), 2, e);
shading interp;
caxis manual
caxis([bottome tope]);
title('Point Level Error');
colorbar;
%nexttile;
subplot(2,2,4);
scatter(X(:,1), X(:,2), 2, ae);
shading interp;
caxis manual
caxis([bottome tope]);
title('Areal Level Error');
colorbar;
set(gcf,'PaperSize',[8 8]);

%saveas(f,"img/ocean.pdf");
print(gcf,'img/ocean.eps','-depsc');
close();
