function [new_qIED, new_qSINGLE, new_cluster_percent] = merge_clusters(NUMBER_OF_CLUSTERS, R, montage_list, el_text, electrodes_with_RAS, clustering_ALL, err)
%% SETTING PARAMS
SEQ = 123;                
SEQ_1 = 123;             
SEQ_2 = 103;             
SEQ_3 = 3;
END = 0;
%% COMBINE INFO from other files and MATCH NEW RAS (mean()) to montage_list
err = flipud(err);
for i = 1:size(err,1)
    del = err(i);
    montage_list(del,:) =[];
    clustering_ALL.qIED(:,del) = [];
    clustering_ALL.qSINGLE(:,del) = [];
    clustering_ALL.qAMPL(:,del) = [];
    clustering_ALL.qAR(:,del) = [];
end

for i = 1:size(montage_list,1)
    
    index_1 = get_index_of_line(i, 1, el_text, montage_list, electrodes_with_RAS);
    index_2 = get_index_of_line(i, 2, el_text, montage_list, electrodes_with_RAS);

    montage_list{i, 3} = get_average('R', electrodes_with_RAS, index_1, index_2);
    montage_list{i, 4} = get_average('A', electrodes_with_RAS, index_1, index_2);
    montage_list{i, 5} = get_average('S', electrodes_with_RAS, index_1, index_2);
end

%% COMPUTE CENTROID 
centroid.qIED = compute_centroid(montage_list, clustering_ALL.qIED, NUMBER_OF_CLUSTERS, clustering_ALL.qAR);
centroid.qSINGLE = compute_centroid(montage_list, clustering_ALL.qSINGLE, NUMBER_OF_CLUSTERS, clustering_ALL.qAR);
centroid.qAMPL = compute_centroid(montage_list, clustering_ALL.qAMPL, NUMBER_OF_CLUSTERS, clustering_ALL.qAR);

%% FILL MERGE_TABLE and R_TABLE

merge_table = fill_merge_table('S', NUMBER_OF_CLUSTERS, R, centroid.qSINGLE);
merge_table = merge_table + fill_merge_table('I', NUMBER_OF_CLUSTERS, R, centroid.qIED);
merge_table = merge_table + fill_merge_table('A', NUMBER_OF_CLUSTERS, R, centroid.qAMPL);

r_table.qIED = fill_r_table(NUMBER_OF_CLUSTERS, R, centroid.qIED);
r_table.qSINGLE = fill_r_table(NUMBER_OF_CLUSTERS, R, centroid.qSINGLE);
r_table.qAMPL = fill_r_table(NUMBER_OF_CLUSTERS, R, centroid.qAMPL);

r_table.average = (r_table.qIED + r_table.qSINGLE + r_table.qAMPL) /3;

%% FIRST RUN - create partly_merged_centroid 

list_of_rep = create_list_of_rep(SEQ, r_table, merge_table);

if isempty(list_of_rep)
    list_of_rep = create_list_of_rep(SEQ_2, r_table, merge_table);
    if isempty(list_of_rep)
        list_of_rep = create_list_of_rep(SEQ_3, r_table, merge_table);
    end
    if isempty(list_of_rep)
        END = 1;
    end
end

if END == 1
    new_qIED = [];
    new_qSINGLE = [];
    new_cluster_percent = [];
else 
    partly_merged_centroid.qIED = create_new_centroid(list_of_rep, centroid.qIED, clustering_ALL);   
    partly_merged_centroid.qSINGLE = create_new_centroid(list_of_rep, centroid.qSINGLE, clustering_ALL);
    partly_merged_centroid.qAMPL = create_new_centroid(list_of_rep, centroid.qAMPL, clustering_ALL);

    % MERGING - result in partly_merged_centroid (sequence: SIA, SI, S) --> GO UNTIL THERES NOTHING TO BE MERGED 
    [partly_merged_centroid.qIED, partly_merged_centroid.qSINGLE, partly_merged_centroid.qAMPL] = go_until_its_possible(R, NUMBER_OF_CLUSTERS, SEQ_1, SEQ_2, SEQ_3, partly_merged_centroid.qSINGLE, partly_merged_centroid.qIED, partly_merged_centroid.qAMPL, clustering_ALL);

    % CREAT q.IED_output cluster_percent_output 
    [qIED_output, qSINGLE_output, cluster_percent_output] = get_result(partly_merged_centroid.qIED, clustering_ALL);
    
    qIED_output_new = [];
    qSINGLE_output_new = [];
    a = 1;
    for col = 1:(size(qIED_output,2)+size(err,1)) %maybe 2 uvidime dal
        if ismember(col,err)
            qIED_output_new = [qIED_output_new nan(size(qIED_output_new,1),1)];
            qSINGLE_output_new = [qSINGLE_output_new nan(size(qSINGLE_output_new,1),1)];
        else
            qIED_output_new = [qIED_output_new qIED_output(:,a)];
            qSINGLE_output_new = [qSINGLE_output_new qSINGLE_output(:,a)];
            a = a + 1;
        end
    end
    
    new_qIED = qIED_output_new;
    new_qSINGLE = qSINGLE_output_new;
    new_cluster_percent = 100*cluster_percent_output/sum(cluster_percent_output);
end
end
function index = get_index_of_line(i, electrode, el_text, montage_list, electrodes_with_RAS)
    part = cell2mat(montage_list(i,2));
    looking_for = el_text{part(electrode)};
    index = find(strcmp(electrodes_with_RAS, looking_for));

end
function average_value = get_average(coord, electrodes_with_RAS, index_1, index_2)
    if coord == 'R'
        coord_index = 2;
    elseif coord == 'A'
        coord_index = 3;
    elseif coord == 'S'
        coord_index = 4;
    end
    
    value_1 = electrodes_with_RAS{index_1, coord_index};
    value_1 = str2double(value_1);
    value_2 = electrodes_with_RAS{index_2, coord_index};
    value_2 = str2double(value_2);
    average_value = mean([value_1 value_2]);
end
function centroid = compute_centroid(montage_list, qNAME, NUMBER_OF_CLUSTERS, qAR)
    centroid = zeros(NUMBER_OF_CLUSTERS, 3);
    
    for j = 1:NUMBER_OF_CLUSTERS
        for k = 1:size(qAR,2) %!!!!!!!!!!!!!!!!!!!montage_list!!!!!!!!!!!!!!!!!!!
            if qAR(j,k) ~= 1
                qNAME(j,k) = 0;
            end
        end
    end
        
    for j = 1:NUMBER_OF_CLUSTERS
        R_sum = 0;
        A_sum = 0;
        S_sum = 0;
        sum_of_qNAME = nansum(qNAME,2);
        for k = 1:size(montage_list,1)
            part_R = cell2mat(montage_list(k,3));
            part_A = cell2mat(montage_list(k,4));
            part_S = cell2mat(montage_list(k,5));
            
            if isnan(qNAME(j,k))
                continue
            else
                R_upper = part_R*qNAME(j,k);
                R_sum = R_sum + R_upper;
                
                A_upper = part_A*qNAME(j,k);
                A_sum = A_sum + A_upper;
                
                S_upper = part_S*qNAME(j,k);
                S_sum = S_sum + S_upper;
            end
        end
        R_centroid = R_sum/sum_of_qNAME(j);
        A_centroid = A_sum/sum_of_qNAME(j);
        S_centroid = S_sum/sum_of_qNAME(j);

        centroid(j,1) = R_centroid;
        centroid(j,2) = A_centroid;
        centroid(j,3) = S_centroid;
    end
    
      %možnost lepšího zápisu
      matice = qNAME ./sum_of_qNAME;
      centroid = matice*cell2mat(montage_list(:,3:5));
end
function X = fill_merge_table(param, num, r, place)
X = zeros(num);
if param == 'S'
     a = 100;
elseif param == 'I'
     a = 20;
elseif param == 'A'
     a = 3;
end

for j = 1:(size(place,1) - 1)
     for k = j+1:size(place,1)
        check = could_be_merged(r, place(k,1), place(k,2), place(k,3), place(j,1), place(j,2), place(j,3));
        if check == 1
            X(j,k) = X(j,k) + a;
        end
     end
    
end
end
function X = fill_merge_table_second(param, num, r, place)
X = zeros(num);
if param == 'S'
     a = 100;
elseif param == 'I'
     a = 20;
elseif param == 'A'
     a = 3;
end

for j = 1:(size(place,1) - 1)
    if cell2mat(place(j,2))==0 && cell2mat(place(j,3))==0 &&  cell2mat(place(j,4))==0
        X(j,:) = -1;
        continue
    end
    for k = j+1:size(place,1)
        if cell2mat(place(k,2))==0 && cell2mat(place(k,3))==0 &&  cell2mat(place(k,4))==0
            X(j,k) = -1;
        else
            check = could_be_merged(r, cell2mat(place(k,2)), cell2mat(place(k,3)), cell2mat(place(k,4)), cell2mat(place(j,2)), cell2mat(place(j,3)), cell2mat(place(j,4)));
            if check == 1
                X(j,k) = X(j,k) + a;
            end
        end
     end
    
end
end
function Y = fill_r_table(num, r, place)
Y = zeros(num);

for j = 1:(size(place,1) - 1)
     for k = j+1:size(place,1)
        [check, r_sqr] = could_be_merged(r, place(k,1), place(k,2), place(k,3), place(j,1), place(j,2), place(j,3));
        if check == 1
            Y(j,k) = r_sqr;
        else
            Y(j,k) = -1; 
        end
     end
    
end
end
function Y = fill_r_table_second(num, r, place)
Y = zeros(num);

for j = 1:(size(place,1) - 1)
     for k = j+1:size(place,1)
        if cell2mat(place(k,2))==0 && cell2mat(place(k,3))==0 &&  cell2mat(place(k,4))==0
            Y(j,k) = -1;
        else
            [check, r_sqr] = could_be_merged(r, cell2mat(place(k,2)), cell2mat(place(k,3)), cell2mat(place(k,4)), cell2mat(place(j,2)), cell2mat(place(j,3)), cell2mat(place(j,4)));
            if check == 1
                 Y(j,k) = r_sqr;
            else
                Y(j,k) = -1; 
            end
        end
     end
    
end
end
function list_of_rep = create_list_of_rep(num, r_table, merge_table)
a = 1;
list_of_rep = [];
for i = 1:size(merge_table,2)
    x = find(merge_table(:,i)>=num);
    if isempty(x)
        continue
    else
        %list_of_rep(1, a) = i;
        for j=1:size(x,1)
            x(j,2) = r_table.average(x(j,1), i);
        end
        value = min(x, [], 1);
        if a >= 2
            %u = find(list_of_rep(2,:)== value(1,1));
             if isempty(find(list_of_rep(2,:)== value(1,1),1)) %u
                 list_of_rep(1, a) = i;
                 list_of_rep(2, a) = value(1,1);
                 a = a + 1; 
             else
                %continue
                break
             end 
         else
             list_of_rep(1, a) = i;
             list_of_rep(2, a) = value(1,1);
             a = a + 1; 
         end
        
    end
end
end
function X = create_new_centroid(list_of_rep, qNAME, clustering_ALL)

for i = 1:size(qNAME,1)
    X{i,1} = i;
    X{i,2} = qNAME(i,1);
    X{i,3} = qNAME(i,2);
    X{i,4} = qNAME(i,3);
end

for i = size(qNAME, 1):-1:1
    w = find(list_of_rep(1,:)== i);
    if isempty(w)
    else
        
        cisla1 = X{list_of_rep(2, w),1};
        cisla2 = X{list_of_rep(1, w),1};
        procento1 = 0;
        procento2 = 0;
        
        for a = 1:size(cisla1,2)
            procento1 = procento1 + clustering_ALL.cluster_percent(cisla1(1,a),1);
        end
        
        for b = 1:size(cisla2,2)
            procento2 = procento2 + clustering_ALL.cluster_percent(cisla2(1,b),1);
        end
        
        
        n = size(X{list_of_rep(2, w),1},2);
        X{list_of_rep(2, w),1} = [X{list_of_rep(2, w),1}, X{i, 1}];
        X{list_of_rep(2, w),2} = (procento1*cell2mat(X(list_of_rep(2, w),2)) + procento2*cell2mat(X(i,2)))/(procento1+procento2);
        X{list_of_rep(2, w),3} = (procento1*cell2mat(X(list_of_rep(2, w),3)) + procento2*cell2mat(X(i,3)))/(procento1+procento2);
        X{list_of_rep(2, w),4} = (procento1*cell2mat(X(list_of_rep(2, w),4)) + procento2*cell2mat(X(i,4)))/(procento1+procento2);
        
        X{i,1} = 0;
        X{i,2} = 0;
        X{i,3} = 0;
        X{i,4} = 0;
    end
        
end
end
function X = create_new_centroid_second(list_of_rep_1, qNAME, clustering_ALL)
X = qNAME;

for i = size(qNAME, 1):-1:1
    w = find(list_of_rep_1(1,:)== i);
    if isempty(w)
    else
        cisla1 = X{list_of_rep_1(2, w),1};
        cisla2 = X{list_of_rep_1(1, w),1};
        procento1 = 0;
        procento2 = 0;
        
        if cisla1 == 0 %kde vznikl error boha? 
            disp('tu');
        end
        
        
        for a = 1:size(cisla1,2)
            procento1 = procento1 + clustering_ALL.cluster_percent(cisla1(1,a),1);
        end
        
        for b = 1:size(cisla2,2)
            procento2 = procento2 + clustering_ALL.cluster_percent(cisla2(1,b),1);
        end
        
        %n = size(X{list_of_rep_1(2, w),1},2);
        X{list_of_rep_1(2, w),1} = [X{list_of_rep_1(2, w),1}, X{i, 1}];
        X{list_of_rep_1(2, w),2} = (procento1*cell2mat(X(list_of_rep_1(2, w),2)) + procento2*cell2mat(X(i,2)))/(procento1+procento2);
        X{list_of_rep_1(2, w),3} = (procento1*cell2mat(X(list_of_rep_1(2, w),3)) + procento2*cell2mat(X(i,3)))/(procento1+procento2);
        X{list_of_rep_1(2, w),4} = (procento1*cell2mat(X(list_of_rep_1(2, w),4)) + procento2*cell2mat(X(i,4)))/(procento1+procento2);
        
        X{i,1} = 0;
        X{i,2} = 0;
        X{i,3} = 0;
        X{i,4} = 0;
    end
        
end
end
function [X, Y, Z] = go_until_its_possible(r, NUMBER_OF_CLUSTERS, seq_1, seq_2, seq_3, qSINGLE, qIED, qAMPL_raw, clustering_ALL)

while(1)
    merge_table_1 = fill_merge_table_second('S', NUMBER_OF_CLUSTERS, r, qSINGLE);
    merge_table_1 = merge_table_1 + fill_merge_table_second('I', NUMBER_OF_CLUSTERS, r, qIED);
    merge_table_1 = merge_table_1 + fill_merge_table_second('A', NUMBER_OF_CLUSTERS, r, qAMPL_raw);

    r_table_1.qIED = fill_r_table_second(NUMBER_OF_CLUSTERS, r, qIED );
    r_table_1.qSINGLE = fill_r_table_second(NUMBER_OF_CLUSTERS, r, qSINGLE);
    r_table_1.qAMPL_raw = fill_r_table_second(NUMBER_OF_CLUSTERS, r, qAMPL_raw);

    r_table_1.average = (r_table_1.qIED + r_table_1.qSINGLE + r_table_1.qAMPL_raw) /3;

    x = find(merge_table_1(:,:)>=seq_1, 1);
    if isempty(x)
        x = find(merge_table_1(:,:)>=seq_2, 1);
        if isempty(x)
            x = find(merge_table_1(:,:)>=seq_3, 1);
            if isempty(x)
                break
            else
                list_of_rep_1 = create_list_of_rep(seq_3, r_table_1, merge_table_1);
            end
        else
            list_of_rep_1 = create_list_of_rep(seq_2, r_table_1, merge_table_1);
        end
    else
        list_of_rep_1 = create_list_of_rep(seq_1, r_table_1, merge_table_1);
    end
    qIED = create_new_centroid_second(list_of_rep_1, qIED, clustering_ALL);
    qSINGLE = create_new_centroid_second(list_of_rep_1, qSINGLE, clustering_ALL);
    qAMPL_raw = create_new_centroid_second(list_of_rep_1, qAMPL_raw, clustering_ALL);
    
    %save('qIED.mat') after every step 
end

    X = qIED;
    Y = qSINGLE;
    Z = qAMPL_raw;

end
function [qIED_output, qSINGLE_output, cluster_percent_output] = get_result(place, clustering)
a = 1;
c = 1;
for i = 1:size(place, 1)
    if cell2mat(place(i,1)) ~= 0
        x = sort(cell2mat(place(i,1)));
        A = zeros(size(x,2),size(clustering.qIED,2));
        C = zeros(size(x,2),size(clustering.qSINGLE,2));
        for j = 1:size(x,2)
            A(j,:) = clustering.qIED(x(j),:);
            C(j,:) = clustering.qSINGLE(x(j),:);
            B(j,1) = clustering.evt_percent(x(j),1); %cluster_percent
        end
    qIED_output(a,:) = sum(A,1);
    qSINGLE_output(c,:) = sum(C,1);
    cluster_percent_output(a,1) = sum(B,1);
    B = [];
    a = a +1;
    c = c +1;
    end
end
end
function [answer,r_sqr] = could_be_merged(r, x, y, z, cx, cy, cz)
    x1 = power((x-cx),2);
    y1 = power((y-cy),2);
    z1 = power((z-cz),2);
    
    r_sqr = x1 + y1 +z1; 
    
    if r_sqr <= power(r,2)
        answer = 1;
    else
        answer = 0; 
    end
end 