%MATLAB code to create Figure 9 in our IEEE Globecom 2008 paper:
%   J.Katto, K.Ogura, Y.Akae, T.Fujikawa, K.Kaneko, S.Zhou: 
%   "Simple Model Analysis and Performance Tuning of Hybrid TCP Congestion 
%   Control", IEEE Globecom 2008, Dec.2008.
%
%   Just type 'GCOM003'. 
%   Default parameter settings are
%      Bandwidth: 10Gbps 
%      RTT: 100msec
%      Packet size: 1500bytes
%   You can modify these parameters and can evaluate several TCP performances.
%
%   Copyright (c) 2008, Katto Lab., Waseda University.
%   All rights reserved.
%
clear all;

% Throughputs & Response Functions

%% Parameter definition
%%% X-axis's Scale
MAX = 11;
%%% Packet size
Ps = 1500 * 8;
%%% Bandwidth
BE = 10000 * 10^6;
%%% RTT
RTT = 100 * 10^(-3);
%%% BDP
W = BE * RTT / Ps;
%%% TCP timer granularity
DMIN = 1 * 10^(-3);
%%% Target number of packets to be stored in a router buffer 
alpha = BE * DMIN / Ps;
%%% Window increase parameter
Winc = BE * DMIN / Ps;
%%% Border probabilities
P1 = 1 / (3/2 * W^2);
P2 = 1 / (3/8 * W^2);
P3 = 1 / (3/2 * (W+alpha)^2);
P4 = 1 / (3/8 * (W+alpha)^2);
P5 = 2*Winc^2 / (3*alpha^2);

%%% Vector Initialization 
scale = 16;
R = zeros(1,MAX*scale);
loss = zeros(1,MAX*scale);
loss2 = zeros(1,MAX*scale);
delay = zeros(1,MAX*scale);
hybrid = zeros(1,MAX*scale);
fusion = zeros(1,MAX*scale);
loss_rtt = zeros(1,MAX*scale);
delay_rtt = zeros(1,MAX*scale);
hybrid_rtt = zeros(1,MAX*scale);
fusion_rtt = zeros(1,MAX*scale);
loss_rtt_ave = zeros(1,MAX*scale);
delay_rtt_ave = zeros(1,MAX*scale);
hybrid_rtt_ave = zeros(1,MAX*scale);
fusion_rtt_ave = zeros(1,MAX*scale);
X = zeros(1,MAX*scale);
b = 1;

%%% Response Functions
delay_response = zeros(1,MAX*scale);
hybrid_response = zeros(1,MAX*scale);
fusion_response = zeros(1,MAX*scale);
reno_response = zeros(1,MAX*scale);
hstcp_response = zeros(1,MAX*scale);
stcp_response = zeros(1,MAX*scale);
reno_response_norm = zeros(1,MAX*scale);

%% Calculate throughputs & response functions
for i=1:MAX*scale
    %%% Parameters
    %%%% PLR (packet loss rate)
    p = 10^(-(MAX-i/scale+1/scale));
    %%%% cwnd
    w = sqrt(8.0/(3.0*p));
    %%%% X-axis
    X(i) = -(MAX-i/scale+1/scale);
    %%%% Timeout penalty
    to = 2 * RTT * (1+p+2*p^2+4*p^3+8*p^4+16*p^5+32*p^6) / (1-p);
    
    %%% Throughputs & response functions
    if (p <= P1) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt(i) = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        delay_num = 1/2 * W * w * Ps;
        delay_rtt(i) = RTT * w/2;
        hybrid_num = 3/8 * w^2 * Ps;
        hybrid_rtt(i) = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        dratio = delay_num / loss_num;
        hratio = hybrid_num / loss_num;
        hybrid_response(i) = log10( (3/4 * w) );
    elseif (p > P1 && p < P2) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt(i) = RTT * w/2 + 1/2 * (w-W)^2 * Ps/BE;
        delay_num = 1/2 * w * W * Ps;
        delay_rtt(i) = RTT * w/2;
        base = 1/2 * (w-W)^2 * Ps;
        hybrid_num = 1/2 * w * W * Ps + base;
        hybrid_rtt(i) = RTT * w/2 + base / BE;
        dratio = delay_num / loss_num;
        hratio = hybrid_num / loss_num;
        hybrid_response(i) = log10( ((w-W)^2/w + W) );
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt(i) = RTT * w/2;
        delay_num = 1/2 * w * W * Ps;
        delay_rtt(i) = RTT * w/2;
        hybrid_num = 1/2 * w * W * Ps;
        hybrid_rtt(i) = RTT * w/2;
        dratio = delay_num / loss_num;
        hratio = hybrid_num / loss_num;
        %hybrid_response(i) = log10( W*(1-p) );
        hybrid_response(i) = log10( hybrid_num / (hybrid_rtt(i) + hratio * to) * RTT/Ps );
    end;
    
    S0 = 1/2 * ((W+alpha)-w/2)^2 * (Winc-1)/Winc;
    S1 = 1/2 * (W-w/2)^2;
    S2 = (W-w/2)*alpha + alpha^2 * (Winc-1) / (2*Winc);
    S3 = 1/2 * (W+alpha-w)^2;
    S4 = 1/2 * alpha^2;
    S5 = 1/2 * (W-w)^2;
    S6 = (W-w) * alpha;
    S7 = 1/8 * (4*W-3*w) * w;
    S8 = 1/8 * w^2 * Winc;
    if ( w > 2*(W+alpha) ) 
        loss2_num = 3/8 * w^2 * Ps;
        loss2_rtt(i) = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        fusion_num = 3/8 * w^2 * Ps;
        fusion_rtt(i) = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        fratio = fusion_num / loss2_num;
        fusion_response(i) = log10( (3/4 * w) );
    elseif ( w > 2*W && w <= 2*(W+alpha) )
        loss2_num = 3/8 * w^2 * Ps;
        loss2_rtt(i) = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        fusion_num = (3/8 * w^2 + S0) * Ps;
        fusion_rtt(i) = RTT * w/2 + (1/8 * (3*w^2 - 4*w*W) + S0) * Ps/BE;
        fratio = fusion_num / loss2_num;
        fusion_response(i) = log10( (3/8 * w^2 + S0) / (w/2) );
    elseif ( w > W+alpha && w <= 2*W )
        loss2_num = 3/8 * w^2 * Ps;
        loss2_rtt(i) = RTT * w/2 + 1/2 * (w-W)^2 * Ps/BE;
        fusion_num = (3/8 * w^2 + S1 + S2) * Ps;
        fusion_rtt(i) = RTT * w/2 + (1/2 * (w-W)^2 + S2) * Ps/BE;
        fratio = fusion_num / loss2_num;
        fusion_response(i) = log10( (3/8 * w^2 + S1 + S2) / (w/2) );
    elseif ( w > W && w <= W+alpha )
        loss2_num = 3/8 * w^2 * Ps;
        loss2_rtt(i) = RTT * w/2 + 1/2 * (w-W)^2 * Ps/BE;
        fusion_num = (3/8 * w^2 + S1 + S2 - S3) * Ps;
        fusion_rtt(i) = RTT * w/2 + (1/2 * (w-W)^2 + S2 - S3) * Ps/BE;
        fratio = fusion_num / loss2_num;
        fusion_response(i) = log10( (3/8 * w^2 + S1 + S2 - S3) / (w/2) );
    elseif ( w > 2*alpha/Winc && w <= W) 
        loss2_num = 3/8 * w^2 * Ps;
        loss2_rtt(i) = RTT * w/2;
        fusion_num = (3/8 * w^2 + S1 - S5 + S2 - S4 - S6) * Ps;
        fusion_rtt(i) = RTT * w/2 + (S2 - S4 - S6) * Ps/BE;
        fratio = fusion_num / loss2_num;
        %fusion_response(i) = log10( (3/8 * w^2 + S1 - S5 + S2 - S4 - S6) / (w/2) );
        fusion_response(i) = log10( fusion_num / (fusion_rtt(i) + fratio * to) * RTT/Ps );
    else
        loss2_num = 3/8 * w^2 * Ps;
        loss2_rtt(i) = RTT * w/2;
        fusion_num = (3/8 * w^2 + S7 + S8) * Ps;
        fusion_rtt(i) = RTT * w/2 + S8 * Ps/BE;
        fratio = fusion_num / loss2_num;
        %fusion_response(i) = log10( (3/8 * w^2 + S7 + S8) / (w/2) );
        fusion_response(i) = log10( fusion_num / (fusion_rtt(i) + fratio * to) * RTT/Ps );
    end

    %%% Throughputs = packets counts / (elapsed time + timeout penalty)
    loss(i) = loss_num / (loss_rtt(i) + to);
    delay(i) = delay_num / (delay_rtt(i) + dratio * to);
    hybrid(i) = hybrid_num / (hybrid_rtt(i) + hratio * to);
    fusion(i) = fusion_num / (fusion_rtt(i) + fratio * to);
    
    loss(i) = loss(i) / 10^6;
    delay(i) = delay(i) / 10^6;
    hybrid(i) = hybrid(i) / 10^6;
    fusion(i) = fusion(i) / 10^6;
    
    %%% Average RTTs
    loss_rtt_ave(i) = loss_rtt(i) * 2/w;
    delay_rtt_ave(i) = delay_rtt(i) * 2/w;
    hybrid_rtt_ave(i) = hybrid_rtt(i) * 2/w;
    fusion_rtt_ave(i) = fusion_rtt(i) * 2/w;

    %%% Response Functions
    delay_response(i) = log10( delay(i) * RTT/Ps * 10^6 );
%    hybrid_response(i) = log10( hybrid(i) * RTT/Ps * 10^6 );
%    fusion_response(i) = log10( fusion(i) * RTT/Ps * 10^6 );
    reno_response(i) = log10( sqrt(1.5/p) );
    reno_response_norm(i) = sqrt(1.5/p) * 1/RTT * Ps / 10^6;
    if p > 10^(-3)
        hstcp_response(i) = reno_response(i);
    else
        hstcp_response(i) = log10( 0.12/p^0.835 );
    end
    stcp_response(i) = log10( 0.08/p );
end;

%% Graphs
%%% PLR-Throughputs
plot(X,loss,'b','LineWidth',2); hold on;
plot(X,delay,'r','LineWidth',2); hold on;
plot(X,hybrid,'g--','LineWidth',2); hold on;
plot(X,fusion,'c--','LineWidth',2); hold on;
%plot(X,reno_response_norm,'b:'); hold on;
line([log10(P1),log10(P1)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P2),log10(P2)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P3),log10(P3)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P4),log10(P4)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P5),log10(P5)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
set(gca,'XLim',[-MAX -1]);
set(gca,'YLim',[0 (BE*1.2)/10^6]);
legend('Loss', 'Delay', 'Hybrid', 'Fusion');
xlabel('Packet loss rate [10\^m]');
ylabel('Throughput [Mbps]');

%%% Response functions
figure;
plot(X,reno_response,'r:','LineWidth',2); hold on;
%plot(X,delay_response,'r','LineWidth',2); hold on;
plot(X,hybrid_response,'g:','LineWidth',2); hold on;
plot(X,fusion_response,'c:','LineWidth',2); hold on;
plot(X,hstcp_response,'b:','LineWidth',2); hold on;
plot(X,stcp_response,'m:','LineWidth',2); hold on;
line([log10(P1),log10(P1)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P2),log10(P2)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P3),log10(P3)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P4),log10(P4)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P5),log10(P5)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
set(gca,'XLim',[-MAX -1]);
set(gca,'YLim',[0 6]);
legend('Reno', 'Hybrid', 'Fusion', 'HSTCP', 'STCP');
xlabel('Packet loss rate [10\^n]]');
ylabel('Packets/RTT [10\^m]');

%%% Average RTTs
figure;
MAX_RTT = RTT*1000*4;
plot(X,loss_rtt_ave*1000,'b','LineWidth',2); hold on;
plot(X,delay_rtt_ave*1000,'r','LineWidth',2); hold on;
plot(X,hybrid_rtt_ave*1000,'g--','LineWidth',2); hold on;
plot(X,fusion_rtt_ave*1000,'c--','LineWidth',2); hold on;
line([log10(P1),log10(P1)],[0,MAX_RTT], 'LineStyle', ':'); 
line([log10(P2),log10(P2)],[0,MAX_RTT], 'LineStyle', ':'); 
line([log10(P3),log10(P3)],[0,MAX_RTT], 'LineStyle', ':'); 
line([log10(P4),log10(P4)],[0,MAX_RTT], 'LineStyle', ':'); 
line([log10(P5),log10(P5)],[0,MAX_RTT], 'LineStyle', ':'); 
set(gca,'XLim',[-MAX -1]);
set(gca,'YLim',[0 MAX_RTT]);
legend('Loss', 'Delay', 'Hybrid', 'Fusion');
xlabel('Packet loss rate [10\^m]');
ylabel('Average RTT [msec]');

%% CSV file outputs
z = zeros(MAX*scale,4);
z(:,1)=X';
z(:,2)=loss';
z(:,3)=delay';
z(:,4)=hybrid';
csvwrite('throughput.txt',z);

z = zeros(MAX*scale,6);
z(:,1)=X';
z(:,2)=reno_response';
z(:,3)=hstcp_response';
z(:,4)=stcp_response';
z(:,5)=hybrid_response';
z(:,6)=fusion_response';
csvwrite('response.txt',z);