clear;clc;close all;
%% Nacitani dat
tabl = readtable('data_ZK.xlsx');
str = 'Jarcová';  % 271
%str = 'Pohořelice'; %301
%str = 'Horní Lideč'; % 365
%str = 'Slavičín'; % 821
%str = 'Holešov'; % 1459
%str = 'Otrokovice'; % 1800
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

% plot(x_val,y_val,'b.')     %Vykresleni vstupnich dat
%% Nastaveni konstant
dmax = 0.15; % 150 metru
cap = 100;
prod = 27/52;
%% Pocitani vzdalenosti
temp = 180/pi;
for i = 1:length(idxs)
    for j = i:length(idxs)
        lat1 = y_val(i)/temp; lat2 = y_val(j)/temp; long1 = x_val(i)/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(i,j) = d;
        D(j,i) = D(i,j);
    end
    b(i) = z_val(i)*prod;
end
D = reshape(D,[],1); D = real(D);
%% Tvorba omezeni
eq21 = sparse(0); eq21(length(idxs),length(idxs)*length(idxs)) = 0;
for i = 1:length(idxs)
    for j = 1:length(idxs)
        eq21(i,(i-1)*length(idxs) + j) = 1;
    end
end
eq21 = sparse([zeros(length(idxs)) eq21]);
eq22 = ones(length(idxs),1);

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

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

f = [ones(1,length(idxs)) zeros(1,length(idxs)*length(idxs))]';
Aineq=[eq31; eq41];
bineq = [eq32; eq42];
Aeq = eq21;
beq = eq22;

%% Vypocet
%% intlinprog
% [x, fval, exitflag, output] = intlinprog (f, 1:(length(D)+length(idxs)), Aineq, bineq, Aeq, beq, zeros(length(D)+length(idxs),1),ones(length(D)+length(idxs),1));
% p = x(1:length(idxs));
% X = x(length(idxs)+1:end); X = reshape(X,length(idxs),length(idxs));

%% cplex
try
    f = [ones(1,length(idxs)) zeros(1,length(idxs)*length(idxs))]';
    Aineq = [eq31; eq41];
    bineq = [eq32; eq42];
    Aeq = eq21;
    beq = eq22;
   
   options = cplexoptimset;
   options.Display = 'on';
   
   [x, fval, exitflag, output] = cplexbilp (f, Aineq, bineq, Aeq, beq, [ ], options);
   p = x(1:length(idxs));
   X = x(length(idxs)+1:end); X = reshape(X,length(idxs),length(idxs));
   fprintf ('\nSolution status = %s\n', output.cplexstatusstring);
   fprintf ('Solution value = %d\n', fval);
catch m
   disp (m.message);
   throw (m);
end
%% Vykresleni
hold on, grid on; axis equal;
plot(x_val, y_val,'k.');
daspect([2 1 1])
for i=1:length(idxs)
    if p(i) > 0.1
        plot(x_val(i), y_val(i),'ro');
        ids = find(X(i,:) == 1); sum(b(ids));
        plot(x_val(ids), y_val(ids),'.')
    end
end
