function [b, E, g, Marg] = bitloadingGA(Hds,Nds,E_total,Gap,PSDmask,const,usedtones, GA_type,bitrate, framerate, Modem);
%% ===========================================================================
% bitloadingGA - Computes bit loading, and energy distribution
%                for a dmt system.
%                         Based on the Levin-Campello algorithm.
% Parameter:    Hds       The channel magnitude response
% Parameter:    Nds	      The noise PSD (mW/Hz)
% Parameter:    E_total	  The total energy (mW) normalized
%                         for subchannel bandwidth
% Parameter:    Gap	      The effective Gap
% Parameter:    PSDmask	  The nominal transmit PSD (mW/Hz)
% Parameter:    const	  The constellation sizes used
% Parameter:    usedtones The dmt tones used
% Parameter:    GA_type	  "0" RA (rate adaptive) / "1" MA ( margin adaptive)
% Parameter:    bitrate	  Bitrate that modem should support for MA algorithm
% Parameter:    framerate Number of frames per second (ADSL=4000, VDSL=??)
% Parameter:    Modem     Modem info (name & to/from nodeno.), for warnings
%                             
%
% Returns:      b    	  Sub-channel bit rate's
% Returns:      E    	  Sub-channel Energy
% Returns:      g         Sub-channel normalized SNR's
% Returns:      Marg      Sub-channel margin's
% Example(s):
%
% Algorithmic details:
%
% Reference:
%  This code is based on LC matlab code from Standford University.
%  Campello J. 'Practical Bit Loading for DMT' IEEE Inter. Conf.
%  on Comm. (ICC), p.796-800.
%  We have modified this algorithm for these constraints:
%    * max. allowed power, max. power for sub-channels (PSDmask),
%    * max. and min. number of bits for subchannel.
%% ===========================================================================

%% ===========================================================================
% Copyright (C):
%       2002-2003 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : FTW's xDSLsimu
% Author(s)     : Driton Statovci (statovci@FTW.at)
%               : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: bitloadingGA.m,v 3.1 2003/12/23 10:35:07 tono Exp $
%% ===========================================================================
% Change History
%      2002-08-12 (DrSt) Created, based on old dmt water-filling code
%      2002-09-24 (DrSt) Fixed initialization of "E_tmp(index)"
%      2003-12-18 (DrSt) Cleaning up  
%% ===========================================================================

Ntotal=length(Hds);                  % Total number of sub-channels

b_min=min(const);                    % Limit constellation sizes
b_max=max(const);

% PSD is multiplied by sqrt(2) (1.5dB) as the PSD is a template
% and the Levin-Campello algorithm expects a peak mask 
PSDmask=PSDmask.*sqrt(2);

% Sub-channel normalized SNR's // Clear out all non-used channels!
g=1e-14*ones(size(Hds));
g(usedtones)=(Hds(usedtones).^2 ./ Nds(usedtones)); 
E=zeros(1,Ntotal);
b=zeros(1,Ntotal);

PSDmask_tmp=PSDmask(usedtones);      % PSD mask of used channels
E_tmp=E(usedtones);           
b_tmp=b(usedtones);

% Used energy so far
E_so_far=0;
% Decision table - QAM signaling
decision_table=Gap./g(usedtones); %
% We suppose that PAM signaling is not used.

if GA_type==0 
    %%%% RA (rate adaptive) bit-loading 
    flag=0;
    while(1)
        [y,index]=min(decision_table);
        if b_tmp(index)==0
            % Increase energy for b_min bit value with the aim to 
            % constraint min. number of bits in sub-channel.
            y= (2^b_min-1)*y; % Increase the used energy with(2^b_min-1)
            flag=1;
        end        
        E_so_far=E_so_far+y;   
        
        if (E_so_far > E_total | sum(b_tmp)>=b_max*length(usedtones));
            E_so_far=E_so_far-y;
            break;
        end
        
        E_tmp(index)=E_tmp(index)+y;
        if (E_tmp(index)>PSDmask_tmp(index)| b_tmp(index)>=b_max)
            E_so_far=E_so_far-y;
            E_tmp(index)=E_tmp(index)-y;
            decision_table(index)=Inf; 
            % No more energy is allowed to be used in this sub-channel
        else   
            if flag==1
                b_tmp(index)=b_tmp(index) + b_min;
                decision_table(index)=(2^b_min)*decision_table(index);
                flag=0;
            else     
                b_tmp(index)=b_tmp(index)+1;
                decision_table(index)=2*decision_table(index);
            end
        end
    end
    % Return to orginal index
    E(usedtones)=E_tmp;
    b(usedtones)=b_tmp;
    
    
else
    %%%% MA (margin adaptive) bit-loading
    flag=0;
    
    bits_per_frame=fix(bitrate*1e3/framerate); % Round to smallest integer
    % This result that the value of calculated bitrate maybe is not the 
    % same as that given as an input parameter. It can be a bit smaller !!!!!

    while(1)
        [y,index]=min(decision_table);
        if b_tmp(index)==0
            y= (2^b_min-1)*y; 
            flag=1;
        end        
        E_so_far=E_so_far+y;   
        
        if (E_so_far > E_total | sum(b_tmp)>= bits_per_frame);
            E_so_far=E_so_far-y;
            break;
        end

        E_tmp(index)=E_tmp(index)+y;    
        if (E_tmp(index)>PSDmask_tmp(index)| b_tmp(index)>=b_max);
            E_so_far=E_so_far-y;
            E_tmp(index)=E_tmp(index)-y;
            decision_table(index)=Inf;
            % No more energy is allowed to be used in this sub-channel
        else   
            if flag==1
                b_tmp(index)=b_tmp(index) + b_min;
                decision_table(index)=(2^b_min)*decision_table(index);
                flag=0;
            else     
                b_tmp(index)=b_tmp(index)+1;
                decision_table(index)=2*decision_table(index);
            end
        end   
    end % End while
    % Return to orginal index
    E(usedtones)=E_tmp;
    b(usedtones)=b_tmp;
        
    % To be calculated only if given bitrate is not supported.
    
    % This part calculates the margin for given bit rate which 
    % is every time negative. 
    if (bits_per_frame>sum(b_tmp))
            
        sprintf('Information for %s modem  LT_Node:%d NT_Node:%d', ...
                Modem.Name, Modem.LT_Node, Modem.NT_Node)
        if (b_max*length(usedtones) < bits_per_frame)
            sprintf('The given bitrate = %d kbps can not be supported for given number of subchannels',bitrate);
            break 
        end
        bitrate_supp=sum(b)*framerate;
        
        E1=zeros(1,Ntotal);
        b1=zeros(1,Ntotal); 
        
        E_tmp=E1(usedtones);
        b_tmp=b1(usedtones);
        
        % Used energy so far
        E_so_far=0;
        % Decision table - QAM signaling
        decision_table=Gap./g(usedtones);
        % We supposed that PAM signaling is not used.
            
        flag=0; % If min. bit in sub-channel is > b_min
        while(1)
            [y,index]=min(decision_table);
            if b_tmp(index)==0
                y= (2^b_min-1)*y; 
                flag=1;
            end
            
            E_so_far=E_so_far+y;
            if (sum(b_tmp)>=bits_per_frame)      
                break;
            end
            if flag==1        
                b_tmp(index)=b_tmp(index) + b_min;
                decision_table(index)=(2^b_min)*decision_table(index);
                flag=0;
            else
                b_tmp(index)=b_tmp(index)+1;
                decision_table(index)=2*decision_table(index);
            end
            E_tmp(index)=E_tmp(index)+y;
        end
        
        Marg = 10*log10(E_total/sum(E_tmp));
        sprintf('To be supported given bitrate = %d kbps this channel need %d dB more Energy',bitrate, -Marg)
        sprintf('For given channel this modem support only the bitrate= %d kbps', bitrate_supp/1e3)
    end % end if bits_per_frame>sum(b_tmp)
end % end LC type

% Calculate the signal margin 
if sum(E)==0
    Marg=0;
else
    Marg = 10*log10(E_total/sum(E));
end