%MATLAB code to create Figure 7(a) 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 'GCOM001'. 
%   Default parameter settings are
%      Bandwidth: 100Mbps
%      RTT: 40msec
%      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;

% Single Flow Model

%% Parameter definition
%%% X-axis's scale
MAX = 6;
%%% Packet size
Ps = 1500 * 8;
%%% Bandwidth
BE = 100 * 10^6;
%%% RTT
RTT = 40 * 10^(-3);
%%% BDP
W = BE * RTT / Ps;
%%% Fusion's Winc
Winc = BE / (10*10^6);
%%% YeAH's STCP weight
yeah_weight = 1.01;
%%% Border probabilities
P1 = 1 / (3/2 * W^2);
P2 = 1 / (3/8 * W^2);

%%% Vector initialization
scale = 16;
loss = zeros(1,MAX*scale);
delay = zeros(1,MAX*scale);
fusion = zeros(1,MAX*scale);
yeah = zeros(1,MAX*scale);
ctcp = zeros(1,MAX*scale);
areno = zeros(1,MAX*scale);
t1 = zeros(1,MAX*scale);
t2 = zeros(1,MAX*scale);
t3 = zeros(1,MAX*scale);
t1_yeah = zeros(1,MAX*scale);
X = zeros(1,MAX*scale);

%% Calculations 
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);
    
    %%% CTCP
    MAXTIME = 10^6;
    %%%% CTCP's time for W (t3)
    w_ctcp2 = zeros(1,MAXTIME);
    w_ctcp2_cwnd = zeros(1,MAXTIME);
    w_ctcp2_dwnd = zeros(1,MAXTIME);
    w_ctcp2(1) = W*0.5;
    w_ctcp2_cwnd(1) = w_ctcp2(1);
    w_ctcp2_dwnd(1) = 0;
    for j=2:MAXTIME
        if ( w_ctcp2(j-1) > W ) 
            w_ctcp2_cwnd(j) = w_ctcp2_cwnd(j-1) * 0.5;
            w_ctcp2_dwnd(j) = w_ctcp2_dwnd(j-1) * 0.5 - w_ctcp2_cwnd(j);
            w_ctcp2(j) = w_ctcp2_dwnd(j) + w_ctcp2_cwnd(j);
            break;
        else
            w_ctcp2_cwnd(j) = w_ctcp2_cwnd(j-1) + 1;
            w_ctcp2_dwnd(j) = w_ctcp2_dwnd(j-1) + 0.125 * w_ctcp2(j-1)^0.75 -1;
            if( w_ctcp2_dwnd(j) < 0) 
                w_ctcp2_dwnd(j) = 0;
            end
            w_ctcp2(j) = w_ctcp2_dwnd(j) + w_ctcp2_cwnd(j);
        end
    end
    index = j;
    t3(i) = MAXTIME;
    for j=2:index
        if( w_ctcp2(j-1) < W && w_ctcp2(j) >= W)
            t3(i) = (j-1) + (W-w_ctcp2(j-1)) / (w_ctcp2(j)-w_ctcp2(j-1)) - 1;
            break;
        elseif( w_ctcp2(1) >= W )
            t3(i) = 0;
            break;
        end
    end
    %%%% CTCP's time for W & w (t1,t2)
    w_ctcp = zeros(1,MAXTIME);
    w_ctcp_cwnd = zeros(1,MAXTIME);
    w_ctcp_dwnd = zeros(1,MAXTIME);
    w_ctcp(1) = w*0.5;
    w_ctcp_cwnd(1) = w_ctcp(1);
    w_ctcp_dwnd(1) = 0;
    w0 = w*0.5;
    for j=2:MAXTIME
        if ( w_ctcp(j-1) > max(W,w) ) 
            w_ctcp_cwnd(j) = w_ctcp_cwnd(j-1) * 0.5;
            w_ctcp_dwnd(j) = w_ctcp_dwnd(j-1) * 0.5 - w_ctcp_cwnd(j);
            w_ctcp(j) = w_ctcp_dwnd(j) + w_ctcp_cwnd(j);
            break;
        else
            w_ctcp_cwnd(j) = w_ctcp_cwnd(j-1) + 1;
            w_ctcp_dwnd(j) = w_ctcp_dwnd(j-1) + 0.125 * w_ctcp(j-1)^0.75 -1;
            if( w_ctcp_dwnd(j) < 0) 
                w_ctcp_dwnd(j) = 0;
            end
            w_ctcp(j) = w_ctcp_dwnd(j) + w_ctcp_cwnd(j);
        end
    end
    index = j;
    t1(i) = MAXTIME;
    t2(i) = MAXTIME;
    for j=2:index
        if( w_ctcp(j-1) < W && w_ctcp(j) >= W)
            t1(i) = (j-1) + (W-w_ctcp(j-1)) / (w_ctcp(j)-w_ctcp(j-1)) - 1;
            break;
        elseif( w_ctcp(1) >= W )
            t1(i) = 0;
            break;
        end
    end
    for j=2:index
        if( w_ctcp(j-1) < w && w_ctcp(j) >= w)
            t2(i) = (j-1) + (w-w_ctcp(j-1)) / (w_ctcp(j)-w_ctcp(j-1)) - 1;
            break;
        elseif( w_ctcp(1) >= w )
            t2(i) = 0;
            break;
        end
    end

    %%% YeAH-TCP
    MAXTIME = 10^6; 
    w_yeah = zeros(1,MAXTIME);
    w_yeah2 = zeros(1,MAXTIME);
    if( W>7*w/8 ) 
        %%%% t1 for YeahTCP
        w_yeah(1) = 7*W/8;
        for j=2:MAXTIME
            if ( w_yeah(j-1) >= W ) 
                break;
            else
                w_yeah(j) = yeah_weight * w_yeah(j-1);
            end
        end
        index = j;
        t1_yeah(i) = MAXTIME;
        for j=2:index
            if( w_yeah(j-1) < W && w_yeah(j) >= W )
                d1 = W - w_yeah(j-1);
                d2 = w_yeah(j) - W;
                t1_yeah(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
    end

    if( W>w )
        %%%% t2 for YeahTCP
        w_yeah2(1) = 7*w/8;
        for j=2:MAXTIME
            if ( w_yeah2(j-1) >= w ) 
                break;
            else
                w_yeah2(j) = yeah_weight * w_yeah2(j-1);
            end
        end
        index = j;
        t2_yeah(i) = MAXTIME;
        for j=2:index
            if( w_yeah2(j-1) < w && w_yeah2(j) >= w )
                d1 = w - w_yeah2(j-1);
                d2 = w_yeah2(j) - w;
                t2_yeah(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
    end

    %%% Parameter Initialization
    loss_num = 0;
    delay_num = 0;
    fusion_num = 0;
    ctcp_num = 0;
    
    %%% Loss(Reno), Delay & Hybrid(Fusion)
    if (W<w0) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        delay_num = 1/2 * W * w * Ps;
        delay_rtt = RTT * w/2;
        fusion_num = 3/8 * w^2 * Ps;
        fusion_rtt = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        dratio = delay_num / loss_num;
        hratio = fusion_num / loss_num;
    elseif (w0<W && w>=W) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/2 * (w-W)^2 * Ps/BE;
        delay_num = 1/2 * w * W * Ps;
        delay_rtt = RTT * w/2;
        base = 1/2 * (w-W)^2 * Ps;
        fusion_num = 1/2 * w * W * Ps + base;
        fusion_rtt = RTT * w/2 + base / BE;
        dratio = delay_num / loss_num;
        hratio = fusion_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
        delay_num = 1/2 * w * W * Ps;
        delay_rtt = RTT * w/2;
        fusion_num = 1/2 * w * W * Ps;
        fusion_rtt = RTT * w/2;
        dratio = delay_num / loss_num;
        hratio = fusion_num / loss_num;
    end;
    loss(i) = loss_num / (loss_rtt + to) / 10^6;
    delay(i) = delay_num / (delay_rtt + dratio * to) / 10^6;
    fusion(i) = fusion_num / (fusion_rtt + hratio * to) / 10^6;

    %%% YeAH-TCP
    if (W<w0) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        yeah_num = 3/8 * w^2 * Ps;
        yeah_rtt = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        yratio = yeah_num / loss_num;
    elseif (w0<W && 7*w/8>=W) 
        base = 1/2 * (w-W)^2 * Ps;
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + base / BE;
        yeah_num = 1/2 * w * W * Ps + base;
        yeah_rtt = RTT * w/2 + base / BE;
        yratio = yeah_num / loss_num;
    elseif (7*w/8<W && W<w) 
        base = 1/2 * (w-W)^2 * Ps;
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + base / BE;
        
        index = fix(t1_yeah(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_yeah(j);
        end
        S1 = S1 + w_yeah(index+1) * (t1_yeah(i)-index);
        S1 = S1 + W*(w/2-t1_yeah(i)) + base / Ps;
        yeah_num = S1 * Ps;
        yeah_rtt = RTT * w/2 + base / BE;
        yratio = yeah_num / loss_num;
    elseif (w<=W && t1_yeah(i) <= w/2)
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;

        index = fix(t1_yeah(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_yeah(j);
        end
        S1 = S1 + w_yeah(index+1) * (t1_yeah(i)-index);
        
        if (S1*Ps < loss_num) 
            yeah_num = loss_num; 
            yeah_rtt = RTT * (t1_yeah(i) + (loss_num/Ps - S1) / W);
            yratio = 1;
        else          
            yeah_num = loss_num;
            yeah_rtt = RTT * t1_yeah(i);
            yratio = 1;
        end        
%       S1 = S1 + W*(w/2-t1_yeah(i));
%       yeah_num = S1 * Ps;
%       yeah_rtt = RTT * w/2;
%       yratio = yeah_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;

        index = fix(t2_yeah(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_yeah2(j);
        end
        S1 = S1 + w_yeah2(index+1) * (t2_yeah(i)-index);
        yeah_num = S1 * Ps;
        yeah_rtt = RTT * t2_yeah(i);
        yratio = yeah_num / loss_num;
    end;
    yeah(i) = yeah_num / (yeah_rtt + yratio * to) / 10^6;
    
    %%% CTCP
    if (W<w0) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        ctcp_num = 1/2 * (w+w0) * (w-w0) * Ps;
        ctcp_rtt = RTT * (w-w0) + (1/2 * (w+w0) * (w-w0) - W*(w-w0)) * Ps/BE;
        cratio = ctcp_num / loss_num;
    elseif (w0<W && w>=W) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/2 * (w-W)^2 * Ps/BE;

        index = fix(t1(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_ctcp(j);
        end
        S1 = S1 + w_ctcp(index+1) * (t1(i)-index);
        S1 = S1 + ((W-w0 - t1(i)) * W) - (1/2*(W+w0)*(W-w0));
        ctcp_num = (1/2*(w+w0)*(w-w0) + S1) * Ps;
        ctcp_rtt = RTT * (w-w0) + (1/2*(w-W)^2) * Ps/BE;
        cratio = ctcp_num / loss_num;
    elseif (w<W && t3(i)<w-w0) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;

        index = fix(t3(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_ctcp(j);
        end
        S1 = S1 + w_ctcp(index+1) * (t3(i)-index);
        
        index = fix(t2(i));
        S2 = 0;
        for j=1:index
            S2 = S2 + w_ctcp(j);
        end
        S2 = S2 + w_ctcp(index+1) * (t2(i)-index);      
%        if (S1*Ps < loss_num) 
%            ctcp_num = loss_num; 
%            ctcp_rtt = RTT * (t3(i) + (loss_num/Ps - S1) / W);
%            cratio = 1;
%        else          
%            ctcp_num = loss_num;
%            ctcp_rtt = RTT * t3(i);
%            cratio = 1;
%        end      
        S1 = S1 + ((w-w0 - t3(i)) * W);
        ctcp_num = S1 * Ps;
        ctcp_rtt = RTT * (w-w0);
        cratio = ctcp_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;

        index = fix(t2(i));
        S2 = 0;
        for j=1:index
            S2 = S2 + w_ctcp(j);
        end
        S2 = S2 + w_ctcp(index+1) * (t2(i)-index);

        ctcp_num = S2 * Ps;
        ctcp_rtt = RTT * t2(i);
        cratio = ctcp_num / loss_num;
    end;    
    ctcp(i) = ctcp_num / (ctcp_rtt + cratio * to) / 10^6;

    %%% ARENO
    if (W<w0) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        areno_num = 3/8 * w^2 * Ps;
        areno_rtt = RTT * w/2 + 1/8 * (3*w^2 - 4*w*W) * Ps/BE;
        aratio = areno_num / loss_num;
    elseif (w0<W && w>=W) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/2 * (w-W)^2 * Ps/BE;
        areno_num = (3/8*w^2 + 1/2*(W+w/2)*(W-w/2)/Winc + W*(W-w/2)*(1-1/Winc) - 1/2*(W+w/2)*(W-w/2)) * Ps;
        areno_rtt = RTT * w/2 + 1/2*(w-W)^2 * Ps/BE;
        aratio = areno_num / loss_num;
    elseif (w<W && W<=Winc*w)
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
        areno_num = (3/8 * W^2/Winc + ((w/2 - W/2/Winc) * W)) * Ps;
        areno_rtt = RTT * w/2;
        aratio = areno_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
        areno_num = 3/8 * W^2/Winc * Ps;
        areno_rtt = RTT * W/2/Winc;
        aratio = areno_num / loss_num;
    end;
    areno(i) = areno_num / (areno_rtt + aratio * to) / 10^6;
end;

%% Graphs
plot(X, loss, X, delay, X, fusion, X, ctcp, X, areno, X, yeah);
line([log10(P1),log10(P1)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P2),log10(P2)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
set(gca,'XLim',[-MAX -1]);
set(gca,'YLim',[0 (BE*1.2)/10^6]);
legend('Loss(Reno)','Delay','Hybrid/Fusion','CTCP','ARENO','YeAH');
xlabel('Packet Loss Rate [10\^n]');
ylabel('Throughput [Mbps]');

%% CSV File Output
z = zeros(5,MAX*scale);
z(1,:)=X;
z(2,:)=loss;
z(3,:)=delay;
z(4,:)=fusion;
z(5,:)=ctcp;
z(6,:)=areno;
z(7,:)=yeah;
csvwrite('single_flow.txt',z');
