function [bridge_out,v_out,deficit_out] = adjust_bridge(bridge_in,v_source,v_target)
% adjust_bridge.m
% 
% bridge   is an  mxn bridging matrix
% v_source is the mx1 source vector
% v_target is the nx1 target vector after bridging

% Initialize variables and display starting info:
bridge = bridge_in;
v_n = bridge'*v_source;
deficit = v_target - v_n;
rms_new = sqrt(mean((deficit).*(deficit))); %root-mean-square deviation represents the sample standard deviation of the residuals
err_start = rms_new;
disp(['RMS error before optimization: ',num2str(err_start)]);
rms_improv = 1;
path_selected = zeros(1,6);

for i=1:200
 

    dummy_deficit = repmat(deficit,1,length(v_target));
    diff_table = dummy_deficit-dummy_deficit';
    [diff_tab_sort,difftabsort_ids] = sort(diff_table(:),'descend');

    max_search =  find(~diff_tab_sort,1)-1;
    direct_transfer   = false;
    indirect_transfer = false;    
    paths = path_selected; % [Col_from,rowlink_from,Col_link,rowlink_to,Row_to,available_transfer]
    for j=1:max_search
        [colA,colB] = ind2sub(size(diff_table),difftabsort_ids(j));
        
        if deficit(colA) > deficit(colB)
            column_TO   = colA;
            column_FROM = colB;
        else
            column_FROM = colA;
            column_TO   = colB;
        end
        
        rowlinks_FROM = find(bridge(:,column_FROM));
        rowlinks_TO = find(bridge(:,column_TO));
        directrowlinks = intersect(rowlinks_FROM,rowlinks_TO);
        if ~isempty(directrowlinks) && deficit(column_FROM)*deficit(column_TO) < 0 %This ensures that the shift is from surplus to deficit
            direct_transfer = true;
            break
        end
            % No direct links found, identify all indirect links (paths)
        for k=1:length(rowlinks_FROM)
            collinks_rowlinks_FROM = find(bridge(rowlinks_FROM(k),:));
            for l=1:length(rowlinks_TO)
                collinks_rowlinks_TO = find(bridge(rowlinks_TO(l),:));
                indirectcollinks = intersect(collinks_rowlinks_FROM,collinks_rowlinks_TO);
                for m=1:length(indirectcollinks)
                    paths(end+1,1:5) = [column_FROM,rowlinks_FROM(k),indirectcollinks(m),rowlinks_TO(l),column_TO];
                end
            end
        end
        if size(paths,1) > 1
            if strcmp(dec2hex(paths(1,1:5)),dec2hex(paths(2,1:5))) && rms_improv<1e-50
                continue
            end
            indirect_transfer = true;
            break
        end
    end

    
    if direct_transfer
        
        transfers = v_source(directrowlinks).*bridge(directrowlinks,column_FROM);
        [max_transfer,max_transfer_id] = max(transfers);
        selected_link = directrowlinks(max_transfer_id);
        transfer_amount = min([max_transfer,deficit(column_TO),-deficit(column_FROM)]);
        transfer_share = transfer_amount/v_source(selected_link);
        bridge(selected_link,column_TO)   = bridge(selected_link,column_TO)   + transfer_share;
        bridge(selected_link,column_FROM) = bridge(selected_link,column_FROM) - transfer_share;
        disp(' ')
        disp(['Iteration ',num2str(i),':'])
        disp(['Transferred ',num2str(transfer_amount),' (',num2str(round(transfer_share*100)),'%) from column ',num2str(column_FROM),' to column ',num2str(column_TO),' via row ',num2str(selected_link)])

    elseif indirect_transfer
        
        if strcmp(dec2hex(paths(1,1:5)),dec2hex(paths(2,1:5))) && rms_improv<1e-50
            continue
            disp(' ')
            disp('No more optimization possible, reduce bridge sparsity to improve match if necessary.')
            i=i-1;
            break
        end

        transfer_required = min(abs(deficit(column_TO)),abs(deficit(column_FROM)));
        for j=2:size(paths,1)
            transfer_available1 = v_source(paths(j,2))*bridge(paths(j,2),paths(j,1));
            transfer_available2 = v_source(paths(j,4))*bridge(paths(j,4),paths(j,3));
            paths(j,6) = min(transfer_available1,transfer_available2);
        end
        
        path_selector = paths(2:end,6);
        path_selected = paths(find(path_selector==max(path_selector),1)+1,:);
        transfer_amount = min(transfer_required,path_selected(6));
        transfer_share1 = transfer_amount/v_source(path_selected(2));
        transfer_share2 = transfer_amount/v_source(path_selected(4));
        
        bridge(path_selected(2),path_selected(1)) = bridge(path_selected(2),path_selected(1)) - transfer_share1;
        bridge(path_selected(2),path_selected(3)) = bridge(path_selected(2),path_selected(3)) + transfer_share1;
        bridge(path_selected(4),path_selected(3)) = bridge(path_selected(4),path_selected(3)) - transfer_share2;
        bridge(path_selected(4),path_selected(5)) = bridge(path_selected(4),path_selected(5)) + transfer_share2;

        bridge(bridge<0) = 0; %Sometimes due to rounding errors the new elements in the bridge can be (very small) negatives, this will cause problems in the optimization.
        
        disp(' ')
        disp(['Iteration ',num2str(i),':'])
        disp(['Transferred ',num2str(transfer_amount),' (',num2str(round(transfer_share1*100)),'%) from column ',num2str(path_selected(1)),' to column ',num2str(path_selected(3)),' via row ',num2str(path_selected(2))])
        disp(['Transferred ',num2str(transfer_amount),' (',num2str(round(transfer_share2*100)),'%) from column ',num2str(path_selected(3)),' to column ',num2str(path_selected(5)),' via row ',num2str(path_selected(4))])

    else
        
        disp(' ')
        disp('No more optimization possible, reduce bridge sparsity to improve match if necessary.')
        i=i-1;
        break
        
    end
    
    v_n = bridge'*v_source;
    deficit = v_target - v_n;
    rms_old = rms_new;
    rms_new =  sqrt(mean(deficit.*deficit));
    rms_improv = rms_old-rms_new;
    disp(['RMS error: ',num2str(rms_new)])

    
end



bridge_out = bridge;
v_out = bridge_out'*v_source;
deficit_out = deficit;
err_stop = rms_new;
disp(['Performed ',num2str(i),' iterations. RMS error was improved from ',num2str(err_start),' to ',num2str(err_stop),'.'])

