clear;
n = 10; %the number of gene trees getting simulated
t1 = 2; 
t2 = 1;
gamma = 0.7;
filename = 'GTs'; %the file which the simulated gene trees will be saved to

fid = fopen(filename,'wt');

for i=1:n
    x1 = exprnd(1);
    if x1 <= t2  %B and C coalesce more recently than hybridization        
        u = rand;  %u decides whether (B,C) goes left or right
        x2 = exprnd(1);
        if x2 <= t1  %(B,C) and A(or D) coalesce more recently than root
            x3 = exprnd(1);
            if u <= gamma  %(B,C) goes left
                gtexp = strcat('(((B:',num2str(x1),',C:',num2str(x1),'):',num2str(t2-x1+x2),',A:',num2str(t2+x2),'):',num2str(t1-x2+x3),',D:',num2str(t1+t2+x3),');');
                senario = 'B and C coalesce within t2, then go left and coalesce with A within t1';
            else
                gtexp = strcat('(((B:',num2str(x1),',C:',num2str(x1),'):',num2str(t2-x1+x2),',D:',num2str(t2+x2),'):',num2str(t1-x2+x3),',A:',num2str(t1+t2+x3),');');
                senario = 'B and C coalesce within t2, then go right and coalesce with D within t1';
            end
        else
            p32 = rand;
            x2 = exprnd(1/3);
            x3 = exprnd(1);
            if p32 <= 1/3  %(B,C) and A coalesce first above root
                gtexp = strcat('(((B:',num2str(x1),',C:',num2str(x1),'):',num2str(t2-x1+t1+x2),',A:',num2str(t2+t1+x2),'):',num2str(x3),',D:',num2str(t2+t1+x2+x3),');');
                senario = 'B and C coalesce within t2, then (B,C) coalesce with A above the root, then with D';
            elseif p32 <= 2/3   %(B,C) and D coalesce first above root
                gtexp = strcat('(((B:',num2str(x1),',C:',num2str(x1),'):',num2str(t2-x1+t1+x2),',D:',num2str(t2+t1+x2),'):',num2str(x3),',A:',num2str(t2+t1+x2+x3),');');
                senario = 'B and C coalesce within t2, then (B,C) coalesce with D above the root, then with A';
            else  %A and D coalesce first above root
                gtexp = strcat('((B:',num2str(x1),',C:',num2str(x1),'):',num2str(t2-x1+t1+x2+x3),',(A:',num2str(t2+t1+x2),',D:',num2str(t2+t1+x2),'):',num2str(x3),');');
                senario = 'B and C coalesce within t2, then coalesce with (A,D) above the root';
            end
        end
    else
        u1 = rand;  %u1 decides whether B goes left or right
        u2 = rand;  %u2 decides whether C goes left or right
        allaboveroot = 0;
        if (u1 <= gamma) && (u2 <= gamma)  %both B and C goes left
            x1 = exprnd(1/3);
            if x1 <= t1
                p32 = rand;
                if p32 <= 1/3  %B and C coalesce first
                    gtexp = strcat('(B:',num2str(t2+x1),',C:',num2str(t2+x1),')');
                    rest = 'A:';
                elseif p32 <= 2/3  %A and C coalesce first
                    gtexp = strcat('(A:',num2str(t2+x1),',C:',num2str(t2+x1),')');
                    rest = 'B:';
                else  %A and B coalesce first
                    gtexp = strcat('(A:',num2str(t2+x1),',B:',num2str(t2+x1),')');
                    rest = 'C:';
                end
                x2 = exprnd(1);
                if x2 <= t1-x1  %A, B and C all coalesce below the root
                    x3 = exprnd(1);
                    gtexp = strcat('((',gtexp,':',num2str(x2),',',rest,num2str(t2+x1+x2),'):',num2str(t1-x1-x2+x3),',D:',num2str(t1+t2+x3),');');
                    senario = 'both B and C goes left, A, B and C all coalesce below the root';
                else
                    p32 = rand;
                    x2 = exprnd(1/3);
                    x3 = exprnd(1);
                    if p32 <= 1/3 %A, B and C coalesce first above the root
                        gtexp = strcat('((',gtexp,':',num2str(t1-x1+x2),',',rest,num2str(t2+t1+x2),'):',num2str(x3),',D:',num2str(t1+t2+x2+x3),');');
                    elseif p32 <= 2/3 % the coalesced node and D coalesce first above the root
                        gtexp = strcat('((',gtexp,':',num2str(t1-x1+x2),',D:',num2str(t2+t1+x2),'):',num2str(x3),',',rest,num2str(t1+t2+x2+x3),');');
                    else
                        gtexp = strcat('(',gtexp,':',num2str(t1-x1+x2+x3),',(D:',num2str(t2+t1+x2),',',rest,num2str(t2+t1+x2),'):',num2str(x3),');');
                    end
                    senario = 'both B and C goes left, two of A, B and C coalesce below the root';
                end
            else
                allaboveroot = 1;
            end
        elseif (u1 > gamma) && (u2 > gamma)  %both B and C goes right
            x1 = exprnd(1/3);
            if x1 <= t1
                p32 = rand;
                if p32 <= 1/3  %B and C coalesce first
                    gtexp = strcat('(B:',num2str(t2+x1),',C:',num2str(t2+x1),')');
                    rest = 'D:';
                elseif p32 <= 2/3  %D and C coalesce first
                    gtexp = strcat('(D:',num2str(t2+x1),',C:',num2str(t2+x1),')');
                    rest = 'B:';
                else  %D and B coalesce first
                    gtexp = strcat('(D:',num2str(t2+x1),',B:',num2str(t2+x1),')');
                    rest = 'C:';
                end
                x2 = exprnd(1);
                if x2 <= t1-x1  %D, B and C all coalesce below the root
                    x3 = exprnd(1);
                    gtexp = strcat('((',gtexp,':',num2str(x2),',',rest,num2str(t2+x1+x2),'):',num2str(t1-x1-x2+x3),',A:',num2str(t1+t2+x3),');');
                    senario = 'both B and C goes right, D, B and C all coalesce below the root';
                else
                    p32 = rand;
                    x2 = exprnd(1/3);
                    x3 = exprnd(1);
                    if p32 <= 1/3 %A, B and C coalesce first above the root
                        gtexp = strcat('((',gtexp,':',num2str(t1-x1+x2),',',rest,num2str(t2+t1+x2),'):',num2str(x3),',A:',num2str(t1+t2+x2+x3),');');
                    elseif p32 <= 2/3 % the coalesced node and A coalesce first above the root
                        gtexp = strcat('((',gtexp,':',num2str(t1-x1+x2),',A:',num2str(t2+t1+x2),'):',num2str(x3),',',rest,num2str(t1+t2+x2+x3),');');
                    else
                        gtexp = strcat('(',gtexp,':',num2str(t1-x1+x2+x3),',(A:',num2str(t2+t1+x2),',',rest,num2str(t2+t1+x2),'):',num2str(x3),');');
                    end
                    senario = 'both B and C goes right, two of D, B and C coalesce below the root';
                end
            else
                allaboveroot = 1;
            end            
        elseif (u1 <= gamma) && (u2 > gamma)
            coal = 0;
            x1 = exprnd(1);
            if x1<=t1 %B and A coalesce below the root
                coal = 1;
            end
            x2 = exprnd(1);
            if x2<=t1 %C and D coalesce below the root
                if coal ==1
                    coal = 3;
                else
                    coal = 2;
                end
            end
            if coal==1
                x2 = exprnd(1/3);
                x3 = exprnd(1);
                p32 = rand;
                if p32<=1/3
                    gtexp = strcat('(((A:',num2str(t2+x1),',B:',num2str(t2+x1),'):',num2str(t1-x1+x2),',C:',num2str(t2+t1+x2),'):',num2str(x3),',D:',num2str(t2+t1+x2+x3),');');
                elseif p32<=2/3
                    gtexp = strcat('(((A:',num2str(t2+x1),',B:',num2str(t2+x1),'):',num2str(t1-x1+x2),',D:',num2str(t2+t1+x2),'):',num2str(x3),',C:',num2str(t2+t1+x2+x3),');');
                else
                    gtexp = strcat('((A:',num2str(t2+x1),',B:',num2str(t2+x1),'):',num2str(t1-x1+x2+x3),',(C:',num2str(t2+t1+x2),',D:',num2str(t2+t1+x2),'):',num2str(x3),');');
                end
                senario = 'B goes left and C goes right; B and A coalesce below the root';
            elseif coal==2
                x1 = x2;
                x2 = exprnd(1/3);
                x3 = exprnd(1);
                p32 = rand;
                if p32<=1/3
                    gtexp = strcat('(((C:',num2str(t2+x1),',D:',num2str(t2+x1),'):',num2str(t1-x1+x2),',A:',num2str(t2+t1+x2),'):',num2str(x3),',B:',num2str(t2+t1+x2+x3),');');
                elseif p32<=2/3
                    gtexp = strcat('(((C:',num2str(t2+x1),',D:',num2str(t2+x1),'):',num2str(t1-x1+x2),',B:',num2str(t2+t1+x2),'):',num2str(x3),',A:',num2str(t2+t1+x2+x3),');');
                else
                    gtexp = strcat('((C:',num2str(t2+x1),',D:',num2str(t2+x1),'):',num2str(t1-x1+x2+x3),',(A:',num2str(t2+t1+x2),',B:',num2str(t2+t1+x2),'):',num2str(x3),');');
                end
                senario = 'B goes left and C goes right; C and D coalesce below the root';
            elseif coal==3
                x3 = exprnd(1);
                gtexp = strcat('((A:',num2str(t2+x1),',B:',num2str(t2+x1),'):',num2str(t1-x1+x3),',(C:',num2str(t2+x2),',D:',num2str(t2+x2),'):',num2str(t1-x2+x3),');');                
                senario = 'B goes left and C goes right; B and A coalesce below the root; C and D coalesce below the root';
            else
                allaboveroot = 1;
            end
        else
            coal = 0;
            x1 = exprnd(1);
            if x1<=t1 %C and A coalesce below the root
                coal = 1;
            end
            x2 = exprnd(1);
            if x2<=t1 %B and D coalesce below the root
                if coal ==1
                    coal = 3;
                else
                    coal = 2;
                end
            end
            if coal==1
                x2 = exprnd(1/3);
                x3 = exprnd(1);
                p32 = rand;
                if p32<=1/3
                    gtexp = strcat('(((A:',num2str(t2+x1),',C:',num2str(t2+x1),'):',num2str(t1-x1+x2),',B:',num2str(t2+t1+x2),'):',num2str(x3),',D:',num2str(t2+t1+x2+x3),');');
                elseif p32<=2/3
                    gtexp = strcat('(((A:',num2str(t2+x1),',C:',num2str(t2+x1),'):',num2str(t1-x1+x2),',D:',num2str(t2+t1+x2),'):',num2str(x3),',B:',num2str(t2+t1+x2+x3),');');
                else
                    gtexp = strcat('((A:',num2str(t2+x1),',C:',num2str(t2+x1),'):',num2str(t1-x1+x2+x3),',(B:',num2str(t2+t1+x2),',D:',num2str(t2+t1+x2),'):',num2str(x3),');');
                end
                senario = 'B goes right and C goes left; B and D coalesce below the root';
            elseif coal==2
                x1 = x2;
                x2 = exprnd(1/3);
                x3 = exprnd(1);
                p32 = rand;
                if p32<=1/3
                    gtexp = strcat('(((B:',num2str(t2+x1),',D:',num2str(t2+x1),'):',num2str(t1-x1+x2),',A:',num2str(t2+t1+x2),'):',num2str(x3),',C:',num2str(t2+t1+x2+x3),');');
                elseif p32<=2/3
                    gtexp = strcat('(((B:',num2str(t2+x1),',D:',num2str(t2+x1),'):',num2str(t1-x1+x2),',C:',num2str(t2+t1+x2),'):',num2str(x3),',A:',num2str(t2+t1+x2+x3),');');
                else
                    gtexp = strcat('((B:',num2str(t2+x1),',D:',num2str(t2+x1),'):',num2str(t1-x1+x2+x3),',(A:',num2str(t2+t1+x2),',C:',num2str(t2+t1+x2),'):',num2str(x3),');');
                end
                senario = 'B goes right and C goes left; B and D coalesce below the root';
            elseif coal==3
                x3 =exprnd(1);
                gtexp = strcat('((A:',num2str(t2+x1),',C:',num2str(t2+x1),'):',num2str(t1-x1+x3),',(B:',num2str(t2+x2),',D:',num2str(t2+x2),'):',num2str(t1-x2+x3),');');                
                senario = 'B goes right and C goes left; C and A coalesce below the root; B and D coalesce below the root';
            else
                allaboveroot = 1;
            end           
        end
        if allaboveroot == 1
            senario = 'all lineages coalesce above the root';
            x1 = exprnd(1/6);
            x2 = exprnd(1/3);          
            x3 = exprnd(1);
            nodes = ['A','B','C','D'];
            randVector = rand(1,4);
            big1 = 1;
            for j=2:4
                if randVector(big1)<randVector(j)
                    big1 = j;
                end
            end
            big2 = 1;
            if big1 == 1
                big2 = 2;
            end
            for j=2:4
                if randVector(big2)<randVector(j) && j~=big1
                    big2 = j;
                end
            end
            newnode = strcat('(',nodes(big1),':',num2str(t1+t2+x1),',',nodes(big2),':',num2str(t1+t2+x1),')');
            rest1 = -1;
            rest2 = -1;
            for j=1:4
                if j~=big1 && j~=big2
                    if rest1==-1
                        rest1 = j;
                    else 
                        rest2 = j;
                    end
                end
            end
            p32 = rand;
            if p32<=1/3
                gtexp = strcat('((',newnode,':',num2str(x2),',',nodes(rest1),':',num2str(t1+t2+x1+x2),'):',num2str(x3),',',nodes(rest2),':',num2str(t1+t2+x1+x2+x3),')');
            elseif p32<=2/3
                gtexp = strcat('((',newnode,':',num2str(x2),',',nodes(rest2),':',num2str(t1+t2+x1+x2),'):',num2str(x3),',',nodes(rest1),':',num2str(t1+t2+x1+x2+x3),')');
            else
                gtexp = strcat('(',newnode,':',num2str(x2+x3),',(',nodes(rest1),':',num2str(t1+t2+x1+x2),',',nodes(rest2),':',num2str(t1+t2+x1+x2),'):',num2str(x3),')');
            end
        end
    end     
    fprintf(fid,'%s\n',gtexp);
end
fclose(fid);
