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

%% População inicial

tamanho_populacao = 30;
numero_geracoes = 100;
populacao_inicial = zeros(tamanho_populacao,M*N);

for t=1:tamanho_populacao
    for m1=1:M
        while 1
            aux1 = randi(N);
            aux2 = randi(K);
            if aux1 >= aux2
                populacao_inicial_aux = randperm(aux1,aux2);
                break
            end
        end
        for j=1:length(populacao_inicial_aux)
            populacao_inicial(t,(m1-1)*N+j) = populacao_inicial_aux(j);
        end
    end
end

%% Opções GA

options = optimoptions('ga',...
    'InitialPopulationMatrix', populacao_inicial,...
    'PopulationSize', tamanho_populacao,...
    'MaxGenerations', numero_geracoes,...
    'EliteCount', 1,...
    'Display', 'off');

%% Função GA

A = []; b = []; Aeq = []; beq = [];
LB = zeros(1,M*N);
UB = zeros(1,M*N) + N;
IntCon=1;
[x,~,~] = ga(@FuncaoObjetivo,M*N,A,b,Aeq,beq,LB,UB,@FuncaoRestricao,IntCon,options);
[~,vazao,~] = FuncaoObjetivo(x);

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

%% Função fitness

    function [vazao_total2,vazao2,vazao_d2d2] = FuncaoObjetivo(gama_celula)
        
        gama_celula2 = reshape(gama_celula.',N,M).';
        vazao2 = zeros(1,M);
        vazao_d2d2 = zeros(1,M);
        for m=1:M
            soma1 = zeros(1,N);
            soma2 = zeros(1,N);
            gama_celula_aux = round(gama_celula2(m,:));
            gama_celula_aux = gama_celula_aux(gama_celula_aux>0);
            for n=1:length(gama_celula_aux)
                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(gama_celula_aux));
            for n=1:length(gama_celula_aux)
                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
            vazao2(1,m) = Wm(1,m)*log2(1 + (Pm(1,m)*Gmb(1,m))/(sigma2_m(1,m) + sum(soma1))) + sum(vazao_n);
            vazao_d2d2(1,m) = sum(vazao_n);
        end
        vazao_total2 = -sum(vazao2./10^6);
    end


%% Restrições

    function [c, ceq] = FuncaoRestricao(gama_celula)
        
        i = 1;
        soma1 = 0;
        soma2 = 0;
        
        gama_celula3 = round(reshape(gama_celula.',N,M).');
        
        for m=1:M
            for n=1:N
                n2 = gama_celula3(m,n);
                if n2>0
                    soma1 = soma1 + Pn(n2) + Gnb(n2);
                    soma2 = soma2 + Pn(n2) + Gnn(n2,n2);
                end
            end
            c(i) = soma1 - Pm(m)*Gmb(m)/gama_m(m) + sigma2_m(m);
            i = i+1;
            for n=1:N
                n2 = gama_celula3(m,n);
                if n2>0
                    c(i) = Pm(m)*Gmn(m,n2) + soma2 - Pn(n2)*Gnn(n2,n2)/gama_n(n2) + sigma2_n(n2);
                    i = i+1;
                end
            end
        end
        ceq = [];
        c = [];
    end

end