% Author: Ian C. Wong and Zukang Shen  
% Copyright (C) 2004   Ian C. Wong
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
% 
% You may reach the author at wongic@mail.utexas.edu.
% Or visit his website at www.ece.utexas.edu/~iwong

% SIPS 2004 Main Simulation File

clear all;
worstpower=1.1565;
N0=worstpower*1e-8;
Ptotal=1;
BER = 1e-3;
Gap = -log(5*BER)/1.6;

channelnum=10;
samplenum=10;
numuser = 8;

epsilong=1e-3;
maxNumOfNewton=0;
maxTotalNewtonIteNum=0;

avenonlintime = zeros(1,numuser);
aveshentime = zeros(1,numuser);
aveiantime = zeros(1,numuser);


diffuservector = [];
iancapamat = [];
shencapamat = [];

cc = 1;

clear shencapa iancapa
N=64;
B=1000000;

noise=B*N0/N;
clear ch


%for different user number
for ii=1:numuser,
    diffuser=2*ii;
    diffuservector(ii)=diffuser;
    K=diffuser;

    totaliancapa=0;
    totalshencapa=0;
    totalnonlincapa = 0;

    shencapavec = zeros(1,diffuser);
    iancapavec = zeros(1,diffuser);

    iannorm = 0;
    shennorm = 0;
    nonlinnorm = 0;


    %for different channel
    for chan=1:channelnum,
        chan
        [env,I,Q]=chtry(K,samplenum,30);
        h = rand(1,K);
        gamma = .064*(h < .5) + .128*((h >= .5) & (h < .8)) + .256*(h>.8);

        %for diff samples
        for diffsamp=1:samplenum,
            diffsamp
            for i=1:K
                user=I(i,:,diffsamp)+sqrt(-1)*Q(i,:,diffsamp);
                ch(i,:)=abs(fft(user,N)).^2/Gap;
            end


            % Shen's subcarrier allocation
            [rheecapa,rheesuballo]=rheesub(Ptotal,ch, N, K, noise, gamma);

            % Shen's power allocation
            t=cputime;  shenp = shenpowerallo(ch,rheesuballo,N,K,Ptotal,noise,gamma); shentime = cputime-t;
            aveshentime(ii) = aveshentime(ii) + shentime;

            % Ian's subcarrier allocation
            [iancapa,iansuballo]=wongsuballo(Ptotal, ch, N, K, noise, gamma);

            % Ian's power allocation
            t=cputime; ianp = wongpowerallo(ch,iansuballo,N,K,Ptotal,noise,gamma); iantime = cputime-t;
            aveiantime(ii) = aveiantime(ii) + iantime;

            for i=1:K,
                shencapa(i) =  waterfilling(shenp(i),rheesuballo(i,:).*ch(i,:)/noise)/N;
                iancapa(i) = waterfilling(ianp(i),iansuballo(i,:).*ch(i,:)/noise)/N;
            end;

            totalshencapa=totalshencapa+sum(shencapa);
            totaliancapa=totaliancapa+sum(iancapa);

            if (chan == 1),
                shencapavec = shencapavec + shencapa;
                iancapavec = iancapavec + iancapa;
                if (chan ==  1 & diffsamp == 1),
                    figure(2);
                    bar([gamma/sum(gamma); iancapavec/sum(iancapavec); shencapavec/sum(shencapavec)]', 'grouped');%; nonlincapavec/sum(nonlincapavec)]', 'grouped');
                    title('snapshot');
                end;
            end;
            iannorm = iannorm + norm(iancapa/sum(iancapa) - gamma/sum(gamma), inf);
            shennorm = shennorm + norm(shencapa/sum(shencapa) - gamma/sum(gamma), inf);
        end

        if (chan == 1),
            iancapavec = iancapavec/(channelnum*samplenum);
            shencapavec = shencapavec/(channelnum*samplenum);
            figure(5);
            bar([gamma/sum(gamma); iancapavec/sum(iancapavec); shencapavec/sum(shencapavec)]', 'grouped');
            legend('Gamma', 'LINEAR', 'ROOT-FINDING');
        end;
        %end diff channel
    end
    maxNumOfNewton;
    maxTotalNewtonIteNum;
    totalshencapavec(ii)=totalshencapa/(channelnum*samplenum);
    totaliancapavec(ii)=totaliancapa/(channelnum*samplenum);
    iannormvec(ii) = iannorm/(channelnum*samplenum);
    shennormvec(ii) = shennorm/(channelnum*samplenum);
end

% total capacities plot
figure(1)
plot(diffuservector,totaliancapavec, 'ko--', diffuservector, totalshencapavec, 'b+-.');%, diffuservector, totalnonlincapavec, 'rx-.');
iansumcapa = sum(totaliancapavec)
shensumcapa = sum(totalshencapavec)

grid on
xlabel('number of users')
ylabel('capacity (bit/s/Hz)')
legend('LINEAR', 'ROOT-FINDING');%, 'NONLIN');
hold off

aveshentime = aveshentime/(channelnum*samplenum)
aveiantime = aveiantime/(channelnum*samplenum)

figure(3);
semilogy(diffuservector,aveiantime, 'ko--', diffuservector,aveshentime, 'b+-.');%, diffuservector, avenonlintime, 'rx-.');
grid on
xlabel('number of users')
ylabel('Ave CPU time (s)')
legend('LINEAR', 'ROOT-FINDING');%, 'NONLIN');
title('Average CPU time comparison');