function plotfigfst()
    global data linewidth1 linewidth2 linewidth3
    close all;
    %ID	L	S	N	K	Py	Beta	Delta	f	Fst	nsuccess	ntotal
    data = sortrows(dlmread('app1.txt','\t'));
    clc;
    
    fig1();%Figure 1. The bias of F̂_IS as a function of the true F_IS
    fig2();%Figure 2. The bias or the RMSE of F ̂_IS as a function of the number L of loci
    %fig3(2);%
    %fig3(4);%
    fig3b(2);%Figure 3. The bias of F̂_IS (or F̂_ST) as a binary function of F_IS and F_ST for Nei’s (1973), Weir & Cockerham’s (1984) and our corrected estimators for diploids. 
    fig3b(4);%Figure S1. The bias of F̂_IS (or F̂_ST) as a binary function of F_IS and F_ST for Nei’s (1973), Weir & Cockerham’s (1984) and our corrected estimators for tetraploids. 
    fig4(0.3);%Figure 4. The bias and the RMSE of F ̂_ST as a function of the number of loci L (F_IS is fixed at 0.3).
    fig5();%Figure 5. The bias of F̂_ST as a function of F_ST
    close all;
end

function fig1()
    w = 12;
    h = 8;
    figure1 = figure('PaperSize',[w h]);
    set(gcf, 'Units','centimeters','Position',[4,1,w,h], 'PaperUnits', 'centimeters', 'PaperPosition', [0,0,w,h])
    
    f1sub(1, 2, 0, 0, figure1);
    f1sub(8, 2, 1, 0, figure1);
    f1sub(1, 4, 0, 1, figure1);
    f1sub(8, 4, 1, 1, figure1);
    
    saveas(gcf, 'f1.eps');
    saveas(gcf, 'f1.pdf');
end

function f1sub(NM, v, x, y, figure1)
    global data 
    
    t1 = data(find(floor(data(:,1) / 1000) == 1 & data(:,2) == v & data(:,4) == 2 & data(:,10) == NM),:);
    t1 = t1(:, [7 14 16 18 30 32 34 20]);
    t1(:,2:8) = t1(:,2:8) - repmat(t1(:,1), [1 7]);
    row = 2;
    col = 2;
    
    top = 0.05;
    bottom = 0.15;
    left = 0.12;
    right = 0.05;
    
    horizontalmargin = 0.012;
    verticalmargin = 0.022;
    
    axes1 = axes('Parent',figure1,...
        'YMinorTick','off',...
        'TickLength',[0.02 0.01],...
        'Position',[left + x * (1 - left - right) / col, bottom + (row - y - 1) * (1 - top - bottom) / row, (1 - left - right - horizontalmargin * col) / col, (1 - top - bottom - verticalmargin * row) / row],...
        'FontName','Times New Roman');
    
    set(gca, 'FontName','Times New Roman');
        
    set(gca, 'YTickLabel',{'', '', '', '', '', '', '', '', ''}, 'YTick', [-0.2:0.05:1]);
    if x == 0
        set(gca, 'YTickLabel',{'-0.2', '', '-0.1', '', '0', '', '0.1', '', '0.2', '', '0.3', '', '4'}, 'YTick', [-0.2:0.05:1]);
    end
    
    set(gca, 'XTickLabel',{'', '', '', '', '', '', '', '', ''}, 'XTick', [0:125:1000])
    if y == row - 1
        set(gca, 'XTickLabel',{'', '', '0.25', '', '0.5', '', '0.75', '', '1'}, 'XTick', [0:0.125:1]);
    end
    
    if y == 0
        if x == 0
            title('Monogamy','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif x == 1
            title('Polygamy','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
    end
    
    if x == 0 && y == 1
        ylabel('\hskip6em Bias($\hat F_{IS}$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
    
    if x == 1 && y == row - 1
        xlabel('$\hskip-9em F_{IS}$', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center');
    end
    
    
    hold on
    t1(:,2) = smooth(t1(:,2), 3);
    t1(:,3) = smooth(t1(:,3), 3);
    t1(:,4) = smooth(t1(:,4), 3);
    t1(:,5) = smooth(t1(:,5), 3);
    t1(:,6) = smooth(t1(:,6), 3);
    t1(:,7) = smooth(t1(:,7), 3);
    t1(:,8) = smooth(t1(:,8), 3);
    
    b = 0.7;
    plot(t1(:,1), t1(:,5), 'k-',  'LineWidth', 2.5, 'Color', [b b b]);
    plot(t1(:,1), t1(:,6), 'k-',  'LineWidth', 2.5, 'Color', [b b b]);
    plot(t1(:,1), t1(:,7), 'k-',  'LineWidth', 2.5, 'Color', [b b b]);
    plot(t1(:,1), t1(:,1) * 0, 'k-',  'LineWidth', 2.5, 'Color', [b b b]);
    plot(t1(:,1), t1(:,2), 'k:',  'LineWidth', 1);
    plot(t1(:,1), t1(:,3), 'k-',  'LineWidth', 1);
    plot(t1(:,1), t1(:,4), 'k-.',  'LineWidth', 1);
    plot(t1(:,1), t1(:,8), 'k--',  'LineWidth', 1);
    
    box on
    Axis = axis;
    Axis(1) = 0;
    Axis(2) = 1;
    Axis(3) = -0.11;
    Axis(4) = .14;
    axis(Axis);
    mid = (Axis(3) + Axis(4))/2;
    
    if x == col - 1
        if y == 0
            ylabel('Diploid','position',[1.1,mid], 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif y == 1
            ylabel('Tetraploid','position', [1.1,mid], 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
    end
end

function fig2()
    w = 16;
    h = 14;
    figure1 = figure('PaperSize',[w h]);
    set(gcf, 'Units','centimeters','Position',[4,1,w,h], 'PaperUnits', 'centimeters', 'PaperPosition', [0,0,w,h])
    
    f2sub(1, 2, 2, 0, 0, figure1);
    f2sub(1, 6, 2, 0, 1, figure1);
    f2sub(8, 2, 2, 0, 2, figure1);
    f2sub(8, 6, 2, 0, 3, figure1);
    
    f2sub(1, 2, 4, 1, 0, figure1);
    f2sub(1, 6, 4, 1, 1, figure1);
    f2sub(8, 2, 4, 1, 2, figure1);
    f2sub(8, 6, 4, 1, 3, figure1);
    
    f2sub(1, 2, 2, 2, 0, figure1);
    f2sub(1, 6, 2, 2, 1, figure1);
    f2sub(8, 2, 2, 2, 2, figure1);
    f2sub(8, 6, 2, 2, 3, figure1);
    
    f2sub(1, 2, 4, 3, 0, figure1);
    f2sub(1, 6, 4, 3, 1, figure1);
    f2sub(8, 2, 4, 3, 2, figure1);
    f2sub(8, 6, 4, 3, 3, figure1);
    
    saveas(gcf, 'f2.eps');
    saveas(gcf, 'f2.pdf');
end

function f2sub(NM, k, v, y, x, figure1)
    global data
    
    t1 = data(find(floor(data(:,1) / 1000) == 2 & data(:,2) == v & data(:,4) == k & data(:,10) == NM),:);
    %bias
    t1b = t1(:, [3 7 14 16 18 42]);
    t1b(:,3:6) = t1b(:,3:6) - repmat(t1b(:,2), [1 4]);
    
    %RMSE
    t1 = t1(:,[3 7 15 17 19 43]);
    t1(:,3:6) = (t1(:,3:6) + t1b(:,3:6) .^ 2) .^ 0.5;
    
    row = 4;
    col = 4;
    
    top = 0.065;
    bottom = 0.088;
    left = 0.088;
    right = 0.025;
    
    horizontalmargin = 0.01;
    verticalmargin = 0.012;
    
    axes1 = axes('Parent',figure1,...
        'YMinorTick','off',...
        'TickLength',[0.02 0.01],...
        'Position',[left + x * (1 - left - right) / col, bottom + (row - y - 1) * (1 - top - bottom) / row, (1 - left - right - horizontalmargin * col) / col, (1 - top - bottom - verticalmargin * row) / row],...
        'FontName','Times New Roman');
    
    set(gca, 'FontName','Times New Roman');
        
    set(gca, 'YTickLabel',{'', '', '', '', '', '', '', '', '', '', '', '', ''}, 'YTick', [0:0.025:0.3]);
    
    if y >= 2 && x == 0
        set(gca, 'YTickLabel',{'0', '', '.05', '', '.1', '', '.15', '', '.2', '', '0.25', '', '.3'}, 'YTick', [0:0.025:0.3]);
    end
    if y < 2 && x == 0
        set(gca, 'YTickLabel',{'-.1', '', '-.05', '', '0'}, 'YTick', [-.1:0.025:0]);
    end
    
    if y == row - 1
        if x == 0 || x == 2
            set(gca, 'XTickLabel',{'', '100', '', '', '250', '', '', '400', '', '', '', '', '', '', ''}, 'XTick', [50:50:500])
        else
            set(gca, 'XTickLabel',{'', '10', '', '', '25', '', '', '40', '', '', '', '', '', ''}, 'XTick', [5:5:50])
        end
    else
        if x == 0 || x == 2
            set(gca, 'XTickLabel',{'', '', '', '', '', '', '', '', '', '', '', '', '', ''}, 'XTick', [50:50:500])
        else
            set(gca, 'XTickLabel',{'', '', '', '', '', '', '', '', '', '', '', '', '', ''}, 'XTick', [5:5:50])
        end
    end
        
    if y == 0
        if x == 0
            title('SNP','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif x == 1
            title('SSR','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif x == 2
            title('SNP','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif x == 3
            title('SSR','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
    end
    
    if x == 0 && y == 1
        ylabel('\hskip6em Bias($\hat F_{IS}$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
    if x == 0 && y == 3
        ylabel('\hskip6.5em RMSE($\hat F_{IS}$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
    
    if x == 1 && y == row - 1
        xlabel('\hskip6em Number of loci ($L$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center');
    end
    
    hold on
    
    len = size(t1, 1) / 3;
    
    for i = 1 : 3
        if y < 2
            tx = t1b([1:len] + (i - 1) * len,:);
        else
            tx = t1([1:len] + (i - 1) * len,:);
        end
        tx(:,3) = smooth(tx(:,3), 5);
        tx(:,4) = smooth(tx(:,4), 5);
        tx(:,5) = smooth(tx(:,5), 5);
        tx(:,6) = smooth(tx(:,6), 5);

        %plot(tx(:,1), tx(:,3), 'k:',  'LineWidth', 0.5);
        color = [0 0 0];
        color(i) = .5;
        plot(tx(:,1), tx(:,4), 'k-',  'LineWidth', 1, 'Color', color);
        plot(tx(:,1), tx(:,5), 'k-.',  'LineWidth', 1, 'Color', color);
        plot(tx(:,1), tx(:,6), 'k--',  'LineWidth', 1, 'Color', color);
    end
    
    box on
    Axis = axis;
    Axis(1) = min(t1(:,1));
    Axis(2) = max(t1(:,1));
    if y >= 2
        Axis(3) = 0;
        Axis(4) = .14;
    else
        Axis(3) = -0.12;
        Axis(4) = 0.02;
    end
    
    axis(Axis);
    xmid = (Axis(1) + Axis(2)) ./ 2;
    mid = (Axis(3) + Axis(4)) ./ 2;
    pos = Axis(1) + (Axis(2) - Axis(1)) * 1.1;
    high = Axis(3) + (Axis(4) - Axis(3)) * 1.23;
    
    if x == col - 1
        if y == 0 || y == 2
            ylabel('Diploid','position',[pos,mid], 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif y == 1 || y == 3
            ylabel('Tetraploid','position', [pos,mid], 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
    end
    
    if y == 0 && x == 0
        text(xmid, high, '\hskip8em Monogamy', 'FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
    if y == 0 && x == 2
        text(xmid, high, '\hskip8em Polygamy', 'FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
end

function fig3(v)
    w = 15;
    h = 9;
    figure1 = figure('PaperSize',[w h]);
    set(gcf, 'Units','centimeters','Position',[4,1,w,h], 'PaperUnits', 'centimeters', 'PaperPosition', [0,0,w,h])
    
    f3sub(1, v, 0, 0, 2, figure1);
    f3sub(1, v, 1, 0, 3, figure1);
    f3sub(1, v, 0, 1, 5, figure1);
    f3sub(1, v, 1, 1, 6, figure1);
    
    saveas(gcf, sprintf('f3_%d.eps', v));
    saveas(gcf, sprintf('f3_%d.pdf', v));
end

function f3sub(NM, v, x, y, id, figure1)
    global data
    
    t1 = data(find(floor(data(:,1) / 1000) == 3 & data(:,2) == v & data(:,4) == 2 & data(:,10) == NM),:);
    t1 = t1(:, [7 8 14 16 18 30 32 34 22 24 26 36 38 40]);
    t1(:,3:8) = t1(:,3:8) - repmat(t1(:,1), [1 6]);
    t1(:,9:14) = t1(:,9:14) - repmat(t1(:,2), [1 6]);
    row = 2;
    col = 2;
    fis = sort(unique(t1(:,1)));
    fst = sort(unique(t1(:,2)));
    
    nx = length(fis);
    ny = length(fst);
    nn = size(t1,1);
    
    tfis1 = zeros(nx, ny);
    tfis2 = zeros(nx, ny);
    tfis3 = zeros(nx, ny);
    tfst1 = zeros(nx, ny);
    tfst2 = zeros(nx, ny);
    tfst3 = zeros(nx, ny);
    
    fis1 = t1(:, [1 2 3]);
    fis2 = t1(:, [1 2 4]);
    fis3 = t1(:, [1 2 5]);
    fst1 = t1(:, [1 2 9]);
    fst2 = t1(:, [1 2 10]);
    fst3 = t1(:, [1 2 11]);
    
    for i = 1 : nn
        xid = find(fis == t1(i,1));
        yid = find(fst == t1(i,2));
        
        tfis1(xid,yid) = t1(i,6);
        tfis2(xid,yid) = t1(i,7);
        tfis3(xid,yid) = t1(i,8);
        tfst1(xid,yid) = t1(i,12);
        tfst2(xid,yid) = t1(i,13);
        tfst3(xid,yid) = t1(i,14);
        
        
         tfis1(xid,yid) = t1(i,3);
         tfis2(xid,yid) = t1(i,4);
         tfis3(xid,yid) = t1(i,5);
         tfst1(xid,yid) = t1(i,9);
         tfst2(xid,yid) = t1(i,10);
         tfst3(xid,yid) = t1(i,11);
    end
    
    top = 0.08;
    bottom = 0.08;
    left = 0.10;
    right = -0.06;
    
    horizontalmargin = 0.07;
    verticalmargin = 0.025;
    
    axes1 = axes('Parent',figure1,...
        'YMinorTick','off',...
        'TickLength',[0.04 0.02],...
        'Position',[left + x * (1 - left - right) / col, bottom + (row - y - 1) * (1 - top - bottom) / row, (1 - left - right - horizontalmargin * col) / col, (1 - top - bottom - verticalmargin * row) / row],...
        'FontName','Times New Roman');
    set(gca, 'FontName','Times New Roman');
    set(gca, 'XTickLabel',{'', '', '.2', '', '.4', '', '.6', '', '.8', '', ''}, 'XTick', [0:0.1:1]);
    set(gca, 'YTickLabel',{'', '', '.2', '', '.4', '', '.6', '', '.8', '', ''}, 'YTick', [0:0.1:1]);
    
    if x == 0
        if y == 0
            zlabel('Bias($\hat F_{IS}$)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        else
            zlabel('Bias($\hat F_{ST}$)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        end
    end
        
    if y == 0
        if x == 0
            title('Nei (1973)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        else
            title('Weir \& Cockerham (1984)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        end
    end
    
    if y >= 1
        xlabel('$F_{ST}$', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        ylabel('$F_{IS}$', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
    end
    %,'HorizontalAlignment','Center','VerticalAlignment','baseline')
     
    hold on
    
    switch id
        case 1, 
            ff = fis1;
            tff = tfis1;
        case 2, 
            ff = fis2;
            tf = tfis2;
            set(gca, 'ZTickLabel',{'-.1', '', '-.06', '', '-.02', '', ''}, 'ZTick', [-.1:.02:0]);
        case 3, 
            ff = fis3;
            tf = tfis3;
            set(gca, 'ZTickLabel',{'-.07', '', '-.04', '', '-.01', '', ''}, 'ZTick', [-.07:.015:-.01]);
        case 4, 
            ff = fst1;
            tf = tfst1;
        case 5, 
            ff = fst2;
            tf = tfst2;
            set(gca, 'ZTickLabel',{'0', '', '.1', '', '.2', '', ''}, 'ZTick', [0:0.05:0.2]);
        case 6, 
            ff = fst3;
            tf = tfst3;
            set(gca, 'ZTickLabel',{'0', '', '.1', '', '.2', '', ''}, 'ZTick', [0:0.05:0.2]);
    end
    
    surf(fis,fst,tf,'Parent',axes1,'EdgeLighting','flat',...
        'FaceLighting','none',...
        'FaceColor','none',...
        'EdgeColor',[0 0 0]);
    %plot3(ff(:,2),ff(:,1),ff(:,3),'k.');
    
    if y == 0
        view(axes1,[-60 50]);
    else
        view(axes1,[-30 50]);
    end
    
    grid(axes1,'on');
    
    Axis = axis;
    Axis(1) = 0;
    Axis(2) = 1;
    Axis(3) = 0;
    Axis(4) = 1;
    axis(Axis);
    zlim(axes1, [min(min(tf)), max(max(tf))]);
end

function fig3b(v)
    w = 14;
    h = 12;
    figure1 = figure('PaperSize',[w h]);
    set(gcf, 'Units','centimeters','Position',[4,1,w,h], 'PaperUnits', 'centimeters', 'PaperPosition', [0,0,w,h])
    
%     f3bsub(1, v, 0, 0, 1, figure1);
    f3bsub(1, v, 0, 0, 2, figure1);
    f3bsub(1, v, 1, 0, 3, figure1);
    f3bsub(1, v, 2, 0, 4, figure1);
    
%     f3bsub(1, v, 0, 1, 11, figure1);
    f3bsub(1, v, 0, 1, 12, figure1);
    f3bsub(1, v, 1, 1, 13, figure1);
    f3bsub(1, v, 2, 1, 14, figure1);
    
    saveas(gcf, sprintf('f3b_%d.eps', v));
    saveas(gcf, sprintf('f3b_%d.pdf', v));
end

function f3bsub(NM, v, y, x, id, figure1)
    global data
    
    t1 = data(find(floor(data(:,1) / 1000) == 3 & data(:,2) == v & data(:,4) == 2 & data(:,10) == NM),:);
    t1 = t1(:, [7 8    14 16 18 20   30 32 34 36    22 24 26 28   36 38 40 42]);
    t1(:,3:10) = t1(:,3:10) - repmat(t1(:,1), [1 8]);
    t1(:,11:18) = t1(:,11:18) - repmat(t1(:,2), [1 8]);
    row = 3;
    col = 2;
    fis = sort(unique(t1(:,1)));
    fst = sort(unique(t1(:,2)));
    
    nx = length(fis);
    ny = length(fst);
    nn = size(t1,1);
    
    tfis1 = zeros(nx, ny);
    tfis2 = zeros(nx, ny);
    tfis3 = zeros(nx, ny);
    tfis4 = zeros(nx, ny);
    tfst1 = zeros(nx, ny);
    tfst2 = zeros(nx, ny);
    tfst3 = zeros(nx, ny);
    tfst4 = zeros(nx, ny);
    
    fis1 = t1(:, [1 2 3]);
    fis2 = t1(:, [1 2 4]);
    fis3 = t1(:, [1 2 5]);
    fis4 = t1(:, [1 2 6]);
    fst1 = t1(:, [1 2 11]);
    fst2 = t1(:, [1 2 12]);
    fst3 = t1(:, [1 2 13]);
    fst4 = t1(:, [1 2 14]);
    
    for i = 1 : nn
        xid = find(fis == t1(i,1));
        yid = find(fst == t1(i,2));
        
        tfis1(xid,yid) = t1(i,7);
        tfis2(xid,yid) = t1(i,8);
        tfis3(xid,yid) = t1(i,9);
        tfis4(xid,yid) = 0;
        
        tfst1(xid,yid) = t1(i,15);
        tfst2(xid,yid) = t1(i,16);
        tfst3(xid,yid) = t1(i,17);
        tfst4(xid,yid) = 0;
    end
    
    top = 0.035;
    bottom = 0.05;
    left = 0.10;
    right = -0.04;
    
    horizontalmargin = 0.07;
    verticalmargin = 0.025;
    
    axes1 = axes('Parent',figure1,...
        'YMinorTick','off',...
        'TickLength',[0.04 0.02],...
        'Position',[left + x * (1 - left - right) / col, bottom + (row - y - 1) * (1 - top - bottom) / row, (1 - left - right - horizontalmargin * col) / col, (1 - top - bottom - verticalmargin * row) / row],...
        'FontName','Times New Roman');
    set(gca, 'FontName','Times New Roman');
    set(gca, 'XTickLabel',{'', '', '.2', '', '.4', '', '.6', '', '.8', '', ''}, 'XTick', [0:0.1:1]);
    set(gca, 'YTickLabel',{'', '', '.2', '', '.4', '', '.6', '', '.8', '', ''}, 'YTick', [0:0.1:1]);
    
    if x == 0
        if y == 4
            zlabel('Wright (1951)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        elseif y == 0
            zlabel('Nei (1973)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        elseif y == 1
            zlabel('Weir \& Cockerham (1984)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        else
            zlabel('Corrected', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        end
    end
        
    if y == 0
        if x == 0
            title('Bias($\hat F_{IS}$)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        else
            title('Bias($\hat F_{ST}$)', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        end
    end
    
    if y == row - 1
        xlabel('$F_{ST}$', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
        ylabel('$F_{IS}$', 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex');
    end
    %,'HorizontalAlignment','Center','VerticalAlignment','baseline')
     
    hold on
    
    switch id
        case 1, 
            ff = fis1;
            tf = tfis1;
            set(gca, 'ZTickLabel',{'0', '', '.2', '', '.4'}, 'ZTick', [0:.1:0.4]);
        case 2, 
            ff = fis2;
            tf = tfis2;
            set(gca, 'ZTickLabel',{'-.1', '', '-.06', '', '-.02', '', ''}, 'ZTick', [-.1:.02:0]);
        case 3, 
            ff = fis3;
            tf = tfis3;
            set(gca, 'ZTickLabel',{'-.07', '', '-.04', '', '-.01', '', ''}, 'ZTick', [-.07:.015:-.01]);
        case 4, 
            ff = fis4;
            tf = tfis4;
            set(gca, 'ZTickLabel',{'-.04', '', '0', '', '.04'}, 'ZTick', [-.04:.02:.04]);
            zlim(axes1, [-.04 0.04]);
        case 11, 
            ff = fst1;
            tf = tfst1;
            set(gca, 'ZTickLabel',{'0', '', '.1', '', '.2', '', ''}, 'ZTick', [0:0.05:0.2]);
        case 12, 
            ff = fst2;
            tf = tfst2;
            set(gca, 'ZTickLabel',{'0', '', '.1', '', '.2', '', ''}, 'ZTick', [0:0.05:0.2]);
        case 13, 
            ff = fst3;
            tf = tfst3;
            set(gca, 'ZTickLabel',{'0', '', '.1', '', '.2', '', ''}, 'ZTick', [0:0.05:0.2]);
        case 14, 
            ff = fst4;
            tf = tfst4;
            set(gca, 'ZTickLabel',{'-.1', '', '0', '', '.1', '', ''}, 'ZTick', [-.1:0.05:0.1]);
            zlim(axes1, [-.1 0.1]);
    end
    
    surf(fis,fst,tf,'Parent',axes1,'EdgeLighting','flat',...
        'FaceLighting','none',...
        'FaceColor','none',...
        'EdgeColor',[0 0 0], 'LineWidth', 0.3);
    plot3(ff(:,2),ff(:,1),ff(:,3),'k.');
    
    if x == 0
        view(axes1,[-60 50]);
    else
        view(axes1,[-30 50]);
    end
    
    grid(axes1,'on');
    
    Axis = axis;
    Axis(1) = 0;
    Axis(2) = 1;
    Axis(3) = 0;
    Axis(4) = 1;
    axis(Axis);
    if y < row - 1
        zlim(axes1, [min(min(tf)), max(max(tf))]);
    end
end

function fig4(fis)
    w = 16;
    h = 14;
    figure1 = figure('PaperSize',[w h]);
    set(gcf, 'Units','centimeters','Position',[4,1,w,h], 'PaperUnits', 'centimeters', 'PaperPosition', [0,0,w,h])
    
    f4sub(fis, 1, 2, 2, 0, 0, figure1);
    f4sub(fis, 1, 6, 2, 0, 1, figure1);
    f4sub(fis, 4, 2, 2, 0, 2, figure1);
    f4sub(fis, 4, 6, 2, 0, 3, figure1);
    f4sub(fis, 1, 2, 4, 1, 0, figure1);
    f4sub(fis, 1, 6, 4, 1, 1, figure1);
    f4sub(fis, 4, 2, 4, 1, 2, figure1);
    f4sub(fis, 4, 6, 4, 1, 3, figure1);
    
    f4sub(fis, 1, 2, 2, 2, 0, figure1);
    f4sub(fis, 1, 6, 2, 2, 1, figure1);
    f4sub(fis, 4, 2, 2, 2, 2, figure1);
    f4sub(fis, 4, 6, 2, 2, 3, figure1);
    f4sub(fis, 1, 2, 4, 3, 0, figure1);
    f4sub(fis, 1, 6, 4, 3, 1, figure1);
    f4sub(fis, 4, 2, 4, 3, 2, figure1);
    f4sub(fis, 4, 6, 4, 3, 3, figure1);
    
    saveas(gcf, sprintf('f4_%f.pdf', fis));
    saveas(gcf, sprintf('f4_%f.eps', fis));
end

function f4sub(fis, NM, k, v, y, x, figure1)
    global data
    
    t1 = data(find(floor(data(:,1) / 1000) == 4 & data(:,2) == v & data(:,4) == k & data(:,7) == fis & data(:,10) == NM),:);
    %bias
    t1b = t1(:, [3 8 22 24 26 44]);
    t1b(:,3:6) = t1b(:,3:6) - repmat(t1b(:,2), [1 4]);
    
    %RMSE
    t1 = t1(:,[3 8 23 25 27 43]);
    t1(:,3:6) = (t1(:,3:6) + t1b(:,3:6) .^ 2) .^ 0.5;
    
    row = 4;
    col = 4;
    
    top = 0.063;
    bottom = 0.088;
    left = 0.08;
    right = 0.025;
    
    horizontalmargin = 0.01;
    verticalmargin = 0.012;
    
    axes1 = axes('Parent',figure1,...
        'YMinorTick','off',...
        'TickLength',[0.02 0.01],...
        'Position',[left + x * (1 - left - right) / col, bottom + (row - y - 1) * (1 - top - bottom) / row, (1 - left - right - horizontalmargin * col) / col, (1 - top - bottom - verticalmargin * row) / row],...
        'FontName','Times New Roman');
    
    set(gca, 'FontName','Times New Roman');
        
    set(gca, 'YTickLabel',{'', '', '', '', '', '', '', '', '', '', '', '', ''}, 'YTick', [0:0.025:0.3]);
    
    if y >= 2 && x == 0
        set(gca, 'YTickLabel',{'0', '', '.05', '', '.1', '', '.15', '', '.2', '', '0.25', '', '.3'}, 'YTick', [0:0.025:0.3]);
    end
    if y < 2 && x == 0
        set(gca, 'YTickLabel',{'0', '', '.05', '', '.1', '', '.15', '', '.2', '', '0.25', '', '.3'}, 'YTick', [0:0.025:0.3]);
    end

    if y == row - 1
        if x == 0 || x == 2
            set(gca, 'XTickLabel',{'', '100', '', '', '250', '', '', '400', '', '', '', '', '', '', ''}, 'XTick', [50:50:500])
        else
            set(gca, 'XTickLabel',{'', '10', '', '', '25', '', '', '40', '', '', '', '', '', ''}, 'XTick', [5:5:50])
        end
    else
        if x == 0 || x == 2
            set(gca, 'XTickLabel',{'', '', '', '', '', '', '', '', '', '', '', '', '', ''}, 'XTick', [50:50:500])
        else
            set(gca, 'XTickLabel',{'', '', '', '', '', '', '', '', '', '', '', '', '', ''}, 'XTick', [5:5:50])
        end
    end
        
    if y == 0
        if x == 0
            title('SNP','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif x == 1
            title('SSR','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif x == 2
            title('SNP','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif x == 3
            title('SSR','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
    end
    
    if x == 0 && y == 1
        ylabel('\hskip6em Bias($\hat F_{ST}$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
    if x == 0 && y == 3
        ylabel('\hskip6.5em RMSE($\hat F_{ST}$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
    
    if x == 1 && y == row - 1
        xlabel('\hskip6em Number of loci ($L$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center');
    end
    
    hold on
    
    len = size(t1, 1) / 3;
    
    for i = 1 : 3
        if y < 2
            tx = t1b([1:len] + (i - 1) * len,:);
        else
            tx = t1([1:len] + (i - 1) * len,:);
        end
        tx(:,3) = smooth(tx(:,3), 5);
        tx(:,4) = smooth(tx(:,4), 5);
        tx(:,5) = smooth(tx(:,5), 5);
        tx(:,6) = smooth(tx(:,6), 5);

        %plot(tx(:,1), tx(:,3), 'k:',  'LineWidth', 0.5);
        color = [0 0 0];
        color(i) = .5;
        %plot(tx(:,1), tx(:,3), 'k:',  'LineWidth', 1, 'Color', color);
        plot(tx(:,1), tx(:,4), 'k-',  'LineWidth', 1, 'Color', color);
        plot(tx(:,1), tx(:,5), 'k-.',  'LineWidth', 1, 'Color', color);
        plot(tx(:,1), tx(:,6), 'k--',  'LineWidth', 1, 'Color', color);
    end
    
    box on
    Axis = axis;
    Axis(1) = min(t1(:,1));
    Axis(2) = max(t1(:,1));
    if y < 2
        Axis(3) = -0.02;
        Axis(4) = 0.12;
    else
        Axis(3) = 0;
        Axis(4) = .12;
    end
    
    axis(Axis);
    xmid = (Axis(1) + Axis(2)) ./ 2;
    mid = (Axis(3) + Axis(4)) ./ 2;
    pos = Axis(1) + (Axis(2) - Axis(1)) * 1.1;
    high = Axis(3) + (Axis(4) - Axis(3)) * 1.23;
    
    if x == col - 1
        if y == 0 || y == 2
            ylabel('Diploid','position',[pos,mid], 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        elseif y == 1 || y == 3
            ylabel('Tetraploid','position', [pos,mid], 'FontSize', 12, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
    end
    
    if y == 0 && x == 0
        text(xmid, high, '\hskip8em Monogamy', 'FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
    if y == 0 && x == 2
        text(xmid, high, '\hskip8em Polygamy', 'FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
    end
end

function fig5()
    global data
    w = 16;
    h = 5;
    figure1 = figure('PaperSize',[w h]);
    set(gcf, 'Units','centimeters','Position',[4,3,w,h], 'PaperUnits', 'centimeters', 'PaperPosition', [0,0,w,h])
    
    row = 1;
    col = 3;
    
    top = 0.01;
    bottom = 0.25;
    left = 0.09;
    right = 0.0;
    
    horizontalmargin = 0.01;
    verticalmargin = 0.012;
    
    x = 0;
    y = 0;
    d1 = data(find(floor(data(:,1) / 1000) == 5 & data(:,2) == 2), [8 24 26 44 46]);
    
    d1(:,2:5) = d1(:,2:5) - repmat(d1(:,1), [1 4]); 
    
    for x = 0:2
        d = d1([x*21+1 : (x+1)*21], :);
        axes1 = axes('Parent',figure1,...
            'YMinorTick','off',...
            'TickLength',[0.02 0.01],...
            'Position',[left + x * (1 - left - right) / col, bottom + (row - y - 1) * (1 - top - bottom) / row, (1 - left - right - horizontalmargin * col) / col, (1 - top - bottom - verticalmargin * row) / row],...
            'FontName','Times New Roman');

        set(gca, 'FontName','Times New Roman');
        plot(axes1, d(:,1), d(:,2), 'k--', 'LineWidth', 1);%Nei
        hold on
        plot(axes1, d(:,1), d(:,3), 'k:', 'LineWidth', 1);%WC
        plot(axes1, d(:,1), d(:,4), 'k-', 'LineWidth', 1);%Huang
        plot(axes1, d(:,1), d(:,5), 'k-.', 'LineWidth', 1);%Hedrick
        box on
        axis([0 1 -0.065 0.14]);
        
        set(gca, 'YTickLabel',{'', '', '', '', '', '', '', '', ''}, 'YTick', [-0.1:0.025:0.15]);
        if x == 0
        set(gca, 'YTickLabel',{'', '', '-.05', '', '0', '', '.05', '', '.1', '', '', '', ''}, 'YTick', [-0.1:0.025:0.15]);
            ylabel('Bias($\hat F_{ST}$)', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
        if x == 1
            xlabel('$F_{ST}$', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center');
        end
        switch x
            case 0
                text(0.9, 0.11, '(A)','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
            case 1
                text(0.9, 0.11, '(B)','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
            case 2
                text(0.9, 0.11, '(C)','FontSize',12,'FontName','Times New Roman', 'Interpreter', 'Latex','HorizontalAlignment','Center','VerticalAlignment','baseline');
        end
    end
    
    saveas(gcf, sprintf('f5.eps'));
    saveas(gcf, sprintf('f5.pdf'));
end
