function [vazao] = AlgoritmoGTMplus(M, K, N, Pm, Pn, gama_m, gama_n, Gnb, Gmb, Gnn, Gmn, sigma2_m, sigma2_n, Wm)

%% Inicialização

Im = (Pm.*Gmb)./gama_m - sigma2_m;
In = (Pn.*Gnb)./gama_n - sigma2_n;

inf_thr = 10;

due_pares = randperm(N, K);
U = [1:M];
gama_celula = cell(1,M);

for n=setdiff(1:N, due_pares)
    m_otimo = WhoGivesMaxUtility(n, U, Pn, Pm, Gnb, Gmb, Gnn, Gmn, Im, Wm, sigma2_m, sigma2_n);
    if m_otimo>0
        gama_celula{m_otimo} = union(gama_celula{m_otimo}, n);
    end
end

idxxxx = 1;
d_set = cell(1,M);
m_dash_vals = zeros(1,M);
while length(U)>0.5
    m_dash = 0;
    leng = -1;

    for i=1:length(U)
        if length(gama_celula{U(i)})>leng
            leng = length(gama_celula{U(i)});
            m_dash = U(i);
        end
    end
    m_dash_vals(idxxxx) = m_dash;

    [cgw, cgn] = FormaGraficoConflito(gama_celula, m_dash, N, Pn, Pm, Gnb, Gmb, Gnn, Gmn, In, Im, Wm, sigma2_m, sigma2_n, inf_thr);

    [temp, mwis] = MaxConjunto(cgw, cgn);
    mwis1 = mwis;

    for i=1:length(mwis)
        n_dash = gama_celula{m_dash}(mwis(i));
        val = Pm(m_dash)*Gmn(m_dash,n_dash) + sum(Pn(setdiff(mwis1,n_dash)).*Gnn(setdiff(mwis1,n_dash), n_dash));
        if val>In(n_dash)
            mwis1 = setdiff(mwis1, mwis(i));
        end
    end
    
    mwis = mwis1;
    for i=1:length(mwis1)
        mwis(i) = gama_celula{m_dash}(mwis1(i));
    end
    [~, idx] = sort(Pn(mwis).*Gmn(m_dash, mwis), 'descend');
    mwis = mwis(idx);
    while sum(Pn(mwis).*Gnb(mwis))>Im(m_dash) && length(mwis)>1.5
        mwis = mwis(2:end);
    end
    
    d_set{m_dash} = mwis;
    temp_set = setdiff(gama_celula{m_dash}, mwis);
    for i=1:length(temp_set)
        m_otimo = WhoGivesMaxUtility(temp_set(i), setdiff(U, m_dash), Pn, Pm, Gnb, Gmb, Gnn, Gmn, Im, Wm, sigma2_m, sigma2_n);
        if m_otimo>0.5
            gama_celula{m_otimo} = union(gama_celula{m_otimo}, temp_set(i));
        end
    end
    idxx = find(U==m_dash);
    U = U([1:(idxx-1) (idxx+1):length(U)]);

    idxxxx = idxxxx + 1;
end

vazao = zeros(1,length(gama_celula));

for m=1:length(gama_celula)
    soma1 = zeros(1,length(cell2mat(gama_celula(1,m))));
    soma2 = zeros(1,length(cell2mat(gama_celula(1,m))));
    if ~isempty(gama_celula(1,m))
        gama_celula_aux = cell2mat(gama_celula(1,m));
        for n=1:length(cell2mat(gama_celula(1,m)))
            soma1(1,n) = Pn(1,gama_celula_aux(n))*Gnb(1,gama_celula_aux(n));
            soma2(1,n) = Pn(1,gama_celula_aux(n))*Gnn(gama_celula_aux(n),gama_celula_aux(n));
        end
        vazao_n = zeros (1,length(cell2mat(gama_celula(1,m))));
        for n=1:length(cell2mat(gama_celula(1,m)))
            vazao_n(1,n) = Wm(1,m)*log2(1 + (Pn(1,gama_celula_aux(n))*Gnn(gama_celula_aux(n),gama_celula_aux(n)))/(sigma2_n(1,gama_celula_aux(n)) + Pm(1,m)*Gmn(m,gama_celula_aux(n)) + sum(soma2)));
        end
    end
    vazao(1,m) = Wm(1,m)*log2(1 + (Pm(1,m)*Gmb(1,m))/(sigma2_m(1,m) + sum(soma1))) + sum(vazao_n);
end

vazao = vazao./10^6; % vazão em mbits/s

end
