clear;clc;close all;
%% Nacitani dat
tabl = readtable('data_ZK.xlsx');
%str = 'Jarcová';  % 271
%str = 'Horní Lideč'; % 365
%str = 'Slavičín'; % 821
str = 'Holešov'; % 1459
%str = 'Otrokovice'; % 1877
%str = 'Kroměříž'; % 4252
idxs = find(strcmp(str,tabl.Vlastn_));

for ind=1:length(idxs)
    y_val(ind) = tabl.Zem_pisn____ka(idxs(ind));
    x_val(ind) = tabl.Zem_pisn_D_lka(idxs(ind));
    z_val(ind) = tabl.Pr_m_rn_Po_etObyvatel(idxs(ind));
end
%% Nastaveni konstant
sz_clust = 300;
dmax = 0.15; % 150 metru
cap = 100;
prod = 27/52;
cl_fval = [];
A = [];
%% Adresy překračující kapacitu
idxs1 = find(z_val > floor(cap/prod)-1);
overload = length(idxs1);
for ind = overload:-1:1
    A{1,ind} = [x_val(idxs1(ind)), y_val(idxs1(ind)), z_val(idxs1(ind))];
    x_val(idxs1(ind)) = []; y_val(idxs1(ind)) = []; z_val(idxs1(ind)) = [];idxs(idxs1(ind)) = [];
    cl_fval(ind) = 1;
    p{1,ind} = 1;
    X{1,ind} = 1;
    b{1,ind} = A{1,ind}(3)*prod';
end
clear idxs1 ind
%% výpočet D
temp = 180/pi;
D = [];
for i = 2:length(idxs)
    for j = i:length(idxs)
        lat1 = y_val(i-1)/temp; lat2 = y_val(j)/temp; long1 = x_val(i-1)/temp; long2 = x_val(j)/temp;
        d = 6371.0*acos((sin(lat1)*sin(lat2))+cos(lat1)*cos(lat2)*cos(long2-long1)); % pocitani opravdove vzdalenosti vzdusnou carou
        D(end+1) = d;
    end
end
clear d;
D = real(D);
%% Clustering
nr_clust = round(length(idxs)/sz_clust);
%Z = linkage(D,'average');
Z = linkage(D,'weighted');
clear D;
%% Vytvořit clustery správné velikosti
max_sz = 2*sz_clust; iter = 0;
while max_sz > sz_clust
    k_idx = cluster(Z,'maxclust',nr_clust + iter);
    sizes = groupcounts(k_idx);
    max_sz = max(sizes);
    iter = iter + 1;
end
clear iter nr_clust;
%% Vytvořit cellarray
for k = 1:length(sizes)
    temp = find(k_idx == k);
    for j = 1:length(temp)
        A_k(j,:) = [x_val(temp(j)),y_val(temp(j)),z_val(temp(j))];
    end
    A{1,overload + k} = A_k;
    A_k = [];
end
clear k;
D = [];
%% Vytvořit rovnice + Cplex
for cluster = overload + 1:length(A)
    %% Vytvořit rovnice
    temp = 180/pi;
    cl_idxs = size(A{1,cluster});
    cl_idxs = cl_idxs(1);
    D = [];
    for i = 1:cl_idxs
        if cl_idxs == 1
            D = 0;
        else
            for j = i:cl_idxs
                lat1 = A{1,cluster}(i,2)/temp; lat2 = A{1,cluster}(j,2)/temp; long1 = A{1,cluster}(i,1)/temp; long2 = A{1,cluster}(j,1)/temp;
                d = 6371.0*acos((sin(lat1)*sin(lat2))+cos(lat1)*cos(lat2)*cos(long2-long1)); % pocitani opravdove vzdalenosti vzdusnou carou
                D(i,j) = d;
                D(j,i) = D(i,j);
            end
        end
        b_k(i) = A{1,cluster}(i,3)*prod;
    end
    D = reshape(D,[],1); D = real(D);
    
    eq21 = []; eq21 = sparse(0); eq21(cl_idxs,cl_idxs*cl_idxs) = 0;
    for i = 1:cl_idxs
        for j = 1:cl_idxs
            eq21(i,(i-1)*cl_idxs + j) = 1;
        end
    end
    eq21 = sparse([zeros(cl_idxs) eq21]);
    eq22 = []; eq22 = ones(cl_idxs,1);

    eq31 = []; eq31 = sparse(-cap*eye(cl_idxs));
    for i = 1:cl_idxs
        eq31 = [eq31 b_k(i)*eye(cl_idxs)];
    end
    eq32 = []; eq32 = 0*ones(cl_idxs,1);

    eq41 = []; eq41 = sparse(zeros(cl_idxs*cl_idxs,cl_idxs));
    eq41_2 = []; eq41_2 = sparse(0); eq41_2(cl_idxs*cl_idxs,cl_idxs*cl_idxs) = 0;
    for i=1:length(D)
        eq41_2(i,i) = D(i);
    end
    eq41 = [eq41, eq41_2];
    eq42 = []; eq42 = dmax*ones(length(D),1);

    f = [ones(1,cl_idxs) zeros(1,cl_idxs*cl_idxs)]';
    Aineq=[eq31; eq41];
    bineq = [eq32; eq42];
    Aeq = eq21;
    beq = eq22;
    %% Výpočet cplex
    try
        options = cplexoptimset;
        options.Display = 'on';
        [x, fval, exitflag, output] = cplexbilp (f, Aineq, bineq, Aeq, beq, [ ], options);
        fprintf ('\nSolution status = %s\n', output.cplexstatusstring);
        fprintf ('Solution value = %d\n', fval);
        p_k = x(1:cl_idxs);
        X_k = x(cl_idxs+1:end); X_k = reshape(X_k,cl_idxs,cl_idxs);
        catch m
            disp (m.message);
            throw (m);
    end
    cl_fval(cluster) = fval;
    p{1,cluster} = p_k;
    X{1,cluster} = X_k;
    b{1,cluster} = b_k';
    b_k = [];
    p_k = [];
    X_k = [];
    
end
full_sol = sum(cl_fval);

fprintf ('Celkove reseni = %d\n', full_sol);
%% Vykresleni
hold on, grid on; axis equal;
plot(x_val, y_val,'k.');
daspect([2 1 1])
for cluster = 1:length(A)
    cl_ind = size(A{1,cluster}); cl_ind = cl_ind(1);
    for i=1:cl_ind
        temp = length(p{1,cluster});
        if temp == 1
            plot(A{1,cluster}(1), A{1,cluster}(2),'ro');
        else
            if p{1,cluster}(i) > 0.1
                plot(A{1,cluster}(i,1), A{1,cluster}(i,2),'ro');
                ids = find(X{1,cluster}(i,:) == 1);
                sum(b{1,cluster}(ids));
                plot(A{1,cluster}(ids,1), A{1,cluster}(ids,2),'.')
            end
        end
    end
end