cd F:/spcluster/revision/supplementary;
rng(22);



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



x = linspace(0.025, 0.975, 50); 
y = linspace(0.025, 0.975, 50);
[xx,yy] = meshgrid(x,y); 
X = [xx(:) yy(:)];




K = 1; % basis size
Theta2 = ones(N, sum(num_basis)+1);
for res = 1:length(num_basis)
    theta = 1/sqrt(num_basis(res))*2.5;
    [knots_s1, knots_s2] = meshgrid(knots_1d{res},knots_1d{res});
    knots = [knots_s1(:), knots_s2(:)];
    for i = 1:num_basis(res)
        disp([res i]);
        d = vecnorm(X -knots(i, :), 2, 2)./theta;
        for j = 1:length(d)
            if d(j) >= 0 & d(j) <= 1
                Theta2(j,i + K) = (1-d(j))^6 * (35 * d(j)^2 + 18 * d(j) + 3)/3;
            else
                Theta2(j,i + K) = 0;
            end
        end
    end
    K = K + num_basis(res);
end


Theta2 = Theta2(:, idx);
Phi2 = Theta2 * Q;

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;

%CC = cov(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  = vecnorm(Eta(i1, :) - Eta(i2, :), 2, 2); %sum(abs(Eta(i1, :) - Eta(i2, :)), 2);  


    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);




    deg = degree(T);
    centers = find(deg >= 4);
    EE = zeros(size(centers,1));
    for i = 1:size(centers,1)
        for j = i:size(centers,1)
             EE(i, j) = findedge(T, centers(i), centers(j));
        end
    end
    rrm = find(sum(EE) > 0);
    centers(rrm) = [];




    dd = pdist2(X(centers, :), X(centers, :));
    dd = dd + 5*eye(size(centers,1));
    CCI = zeros(size(centers,1)-1, 2);
    for i= 1: size(centers,1)-1
        [~, ij] = min(dd(i,i:end));
        ij = i - 1 + ij;
        sp = shortestpath(T, centers(i), centers(ij));
        if sum(ismember(sp(2:end-1), centers))
            if ij <= size(centers,1)
                ij1 = min(find(ismember(sp(2:end-1), centers)));
                ij = find(centers == sp(ij1+1));
            end
        end
        CCI(i, :) = [ centers(i), centers(ij) ];
    end


    IDX = [];
    for i= 1: size(CCI,1)-1 
        U = shortestpathtree(T, CCI(i,1), CCI(i,2));
        [~, maxidx] = max(U.Edges.Weight);
        IDX = [IDX findedge(T, U.Edges.EndNodes(maxidx, 1), U.Edges.EndNodes(maxidx, 2))];
    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) );
 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  = vecnorm(Eta(i1, :) - Eta(i2, :), 2, 2); %sum(abs(Eta(i1, :) - Eta(i2, :)), 2);  


    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);




    deg = degree(T);
    centers = find(deg >= 4);
    EE = zeros(size(centers,1));
    for i = 1:size(centers,1)
        for j = i:size(centers,1)
             EE(i, j) = findedge(T, centers(i), centers(j));
        end
    end
    rrm = find(sum(EE) > 0);
    centers(rrm) = [];




    dd = pdist2(X(centers, :), X(centers, :));
    dd = dd + 5*eye(size(centers,1));
    CCI = zeros(size(centers,1)-1, 2);
    for i= 1: size(centers,1)-1
        [~, ij] = min(dd(i,i:end));
        ij = i - 1 + ij;
        sp = shortestpath(T, centers(i), centers(ij));
        if sum(ismember(sp(2:end-1), centers))
            if ij <= size(centers,1)
                ij1 = min(find(ismember(sp(2:end-1), centers)));
                ij = find(centers == sp(ij1+1));
            end
        end
        CCI(i, :) = [ centers(i), centers(ij) ];
    end


    IDX = [];
    for i= 1: size(CCI,1)-1 
        U = shortestpathtree(T, CCI(i,1), CCI(i,2));
        [~, maxidx] = max(U.Edges.Weight);
        IDX = [IDX findedge(T, U.Edges.EndNodes(maxidx, 1), U.Edges.EndNodes(maxidx, 2))];
    end
    
    
    
U = T.rmedge(IDX);
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, :), "all");
end

e  = abs(Yhat - mean(Y2,2));
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([Y2, Yhat], [], 'all');
topY    = max([Y2, Yhat], [], 'all');

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



f = tiledlayout(2,2);
nexttile;
scatter(XX(:,1), XX(:,2), 1, YY);
shading interp;
caxis manual
caxis([bottomY topY]);
title('Point Level Data');
colorbar;
nexttile;
scatter(X(:,1), X(:,2), 1, Yhat);
shading interp;
caxis manual
caxis([bottomY topY]);
title('Areal Level Data');
colorbar;
nexttile;
scatter(X(:,1), X(:,2), 1, e);
shading interp;
caxis manual
caxis([bottome tope]);
title('Point Level Error');
colorbar;
nexttile;
scatter(X(:,1), X(:,2), 1, ae);
shading interp;
caxis manual
caxis([bottome tope]);
title('Areal Level Error');
colorbar;



print(gcf,'img/sim.eps','-depsc');
close();

