% letter chi 4

% cagingfig show different metrics

close all;
linewidth = 1.25;
smoothflag = true;
% set up figure
widthpix = 375;
heightpix = 300;

leftmargin = 67;
bottommargin = 80;
topmargin = 85;
rightmargin = 90;
centermargin = 0;
widthtotal = widthpix*2+leftmargin+rightmargin;
heighttotal = heightpix*2+topmargin+bottommargin;
left = leftmargin/widthtotal;
bottom = bottommargin/heighttotal;
center = centermargin/widthtotal;
width = widthpix/widthtotal;
height = heightpix/heighttotal;

psize = 26; % font size
msize = 9; % marker size
fig = figure('rend','painters','pos',[100 100 widthtotal heighttotal]);

%% panel (a) chi with time
panelA = axes('Position', [left, bottom+height, width, height]);

phihat = 5;
dims = {'4D', '8D', '12D', '16D', '20D'};

n=length(dims);
rpp=(1:n)';
colors = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];
for rp=1:length(dims)
    plot(nan, nan, 'color', colors(rp,:), 'linewidth', 1.5);
    if 1==rp, hold on, end
end

for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    ndimsqr = ndim*ndim*ndim;
    phistr = num2str(phihat*ndim);
    fname = strcat('../data/rlgdynamics/', dim, '/', phistr, '/msd.dat');
    data = dlmread(fname, '', 1, 0);
    
    xs = data(:,1)*sqrt(ndim);
    ys = (data(:,3) ./ (data(:,2).*data(:,2)) - 1 )*ndim;
    ys2 = ( (data(:,3)-data(:,4) ) ./ (data(:,2).*data(:,2))  )*ndim;
    plot(xs, ys, 'color', colors(rp,:) );
    plot(xs, ys2, '--', 'color', colors(rp,:) );
    if phihat < 5
        meanys = mean(ys(floor(end/2):end));
    else
        meanys = mean(ys(end-4:end));
    end
    %fprintf('%g\t%g\t%g\n', phihat, ndim, meanys);
end

set(gca, 'xscale', 'log');
xlabh = xlabel('$\hat{t}$', 'interpreter', 'latex');
set(xlabh, 'Units', 'Normalized', 'Position', [0.5,1.13, 0]);
ylabh = ylabel('$\hat\chi^\mathrm{S}/\hat\Delta^2$', 'interpreter', 'latex');
set(ylabh, 'Units', 'Normalized', 'Position', [-0.07, 0.5, 0]);

h = legend(dims, 'interpreter', 'latex', 'location', 'southeast', 'box', 'off');
set(h,'Position', get(h, 'Position')+ [-0.01 0.048 0 0]);

xlim([1e-1,1e5]);
ylim([0, 10]);
xticks([1e0,1e2,1e4]);
yticks([0, 5, 10]);
set(gca, 'XAxisLocation', 'top');

%legend(dims, 'location', 'northwest');
%legend boxoff;
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

box on;

%% panel (b) chi_th distribution
panelB = axes('Position', [left+width, bottom+height, width, height]);


distributionfile = 'Dtdist_log.dat';
fitdata = dlmread('../data/rlgdynamics/chis.dat', '', 1, 0);
holdflag = false;
smoothflag = true;

for rp = 1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    phi = ndim*phihat;
    phistr = num2str(phi, '%.6g');
    foldername = strcat('../data/rlgdynamics/', dim, '/', phistr, '/');
    buffname = strcat(foldername, distributionfile);
    M = dlmread(buffname);
    % deal with outliers
    idxs = find(M(:,2)==0);
    [a_, imax] = max( diff(idxs) );
    xs = M(:,1);
    ys = M(:,2);
    ys(1:idxs(imax)) = 0;
    ys(idxs(imax+1):end) = 0;
    dx = log(xs(2)) - log(xs(1));
    edges = log(xs(1)) - dx/2 + ( 0:numel(xs) ).'*dx;
    if smoothflag % only smooth the tail
         ys = smooth(ys);
    end
    % normalization
    ys = ys ./ ( diff( exp(edges) ) * sum(ys) );
    plot(xs, ys, 'color', colors(rp,:) );
    if ~holdflag
        holdflag = true;
        hold on;
    end
    % fit line
    fitline = fitdata(fitdata(:,1)==ndim & fitdata(:,2)==phi, :);
    ys = exp(-(log(xs/ndim) - fitline(6) ).^2 ./ (2*fitline(7)*fitline(7))) ./ (xs *fitline(7) * sqrt(2*pi));
    plot(xs, ys, ':', 'color', colors(rp,:));
end
mfvals = dlmread('../data/Deltamf.dat');
yval = interp1(mfvals(:,2), mfvals(:, 1), phihat, 'spline');

plot([yval, yval], [1e-3, 1e1], '--k', 'linewidth', 0.75);
set(gca, 'xscale', 'log', 'yscale', 'log');
xlim([4e-3, 3]);
ylim([1e-3, 1e1]);
xticks([1e-2, 1e-1, 1e0]);
yticks([1e-3, 1e-1, 1e1]);
xlabel('$\hat{\Delta}$', 'interpreter', 'latex');
ylabel('$P(\hat{\Delta})$', 'interpreter', 'latex');

set(gca, 'XAxisLocation', 'top', 'YAxisLocation', 'right');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (c) chi_th vs cage size

panelC = axes('Position', [left, bottom, width, height]);

dims = {'4D', '8D', '12D', '16D', '20D'};
phihats = [3, 3.5, 4, 5, 7, 10, 15, 20, 30, 40];
holdflg = true;

n=length(dims);
rpp=(1:n)';
colors = [(rpp-1)/(n-1), zeros(n,1), 1 - (rpp-1)/(n-1)];

for rp=1:length(dims)
    plot(nan, nan, 'color', colors(rp,:) );
    if 1==rp, hold on, end
end
lgds = dims;

mftchith = dlmread('../data/MFTChi/Chith.dat');  % is chi_th
mftchhet = dlmread('../data/MFTChi/Chihet.dat'); % is chi_het
allallresults = {};
for rp=1:length(dims)
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    
    % filter out outliers
    thisdim = fitdata(:,1)==ndim;
    csize = fitdata(thisdim, 3); chith=fitdata(thisdim, 4); chihet2=fitdata(thisdim, 5);
    csizesqr = csize.*csize; 
    allresults = [fitdata(thisdim,2)/ndim, csize*ndim, ndim*(chith-chihet2)./csizesqr, ndim*chihet2./csizesqr ];

    plot(allresults(:,2), allresults(:, 3), 's-', 'markersize', msize, 'color', colors(rp,:) );
    allallresults{end+1} = allresults;
    if 1==rp, hold on, end
end

% plot theoretical line
plot(mftchith(:,1), mftchith(:,2), '-', 'color', 'k');
ylabh = ylabel('$\hat\chi^\mathrm{S}_\mathrm{th} / \hat\Delta^2$', 'interpreter', 'latex');
set(ylabh, 'Units', 'Normalized', 'Position', [-0.06, 0.5, 0]);

xlim([1e-3, 2e0]);
ylim([2.6, 23]);
xticks([1e-3, 1e-2, 1e-1]);
yticks([3, 5, 10 ]);
xlabel('$\hat\Delta$', 'interpreter', 'latex');

set(gca, 'xscale', 'log', 'yscale', 'log');
%set(gca, 'yscale', 'log');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);

%% panel (c) inset chi_th correction

insetC = axes('Position', [left+0.05, bottom+0.2, width*0.48, height*0.47]);

for rp=1:length(dims)
    ndim = str2double(dims{rp}(1:end-1));
    allresults = allallresults{rp};
    xs = allresults(:,2);
    ybase = interp1(mftchith(:,1), mftchith(:,2), allresults(:,2));
    ys = (allresults(:, 3) - ybase)./ybase;
    plot(xs, ndim*ys, 's-', 'markersize', msize-1, 'color', colors(rp,:) );
    if 1==rp, hold on, end
end

xlim([1e-3, 1e0]);
ylim([6e-1, 2e1]);
xticks([1e-3, 1e-2, 1e-1, 1e0]);
yticks([1, 4, 10]);
set(gca, 'xscale', 'log', 'yscale', 'log');
%xlabel('$\hat\Delta$', 'interpreter', 'latex');
ylabel('$d \times \delta \hat\chi^\mathrm{S}_\mathrm{th}$', 'interpreter', 'latex');
set(gca, 'YAxisLocation', 'right');

set(gca, 'fontname', 'times new roman', 'fontsize', psize-2);

%% panel (d) chi_het


panelD = axes('Position', [left+width, bottom, width, height]);

for rp=1:length(dims)
        plot(nan, nan, 'color', colors(rp,:) );
        if 1==rp, hold on, end
end


for rp=[1,2,4]
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    fts = fitdata(:,1) == ndim;
    phihats = fitdata(fts, 2)/ndim;
    csize = exp(fitdata(fts, 6));
    typcsize = exp(fitdata(fts, 6) + fitdata(fts, 7).^2/2);
    typhet = (exp(fitdata(fts, 7).^2) - 1 ) .* exp(2*fitdata(fts, 6) + fitdata(fts, 7).^2);
    plot(csize*ndim, ndim*typhet ./ (typcsize.*typcsize), ':x', 'color', colors(rp,:));
    if 1==rp, hold on, end
end
 

for rp=1:length(dims)
    allresults = allallresults{rp};
    xs = allresults(:,2);
    ybase = interp1(mftchith(:,1), mftchith(:,2), allresults(:,2));
    ys = (allresults(:, 3) - ybase)./ybase;
    plot(allresults(:,2), allresults(:, 4), 'x-', 'markersize', msize+1, 'color', colors(rp,:) );
end


plot(mftchhet(:,1), mftchhet(:,2), '--', 'color', 'k' );

ylabh=ylabel('$\hat\chi^\mathrm{S}_\mathrm{het} / \hat\Delta^2$', 'interpreter', 'latex');
set(ylabh, 'Units', 'Normalized', 'Position', [1.14, 0.5, 0]);
xlim([1e-3, 2e0]);
ylim([2.6, 23]);
xticks([1e-3, 1e-2, 1e-1, 1e0]);
yticks([3, 5, 10]);
xlabel('$\hat\Delta$', 'interpreter', 'latex');

set(gca, 'xscale', 'log', 'yscale', 'log');
set(gca, 'yAxisLocation', 'right');
set(gca, 'fontname', 'times new roman', 'fontsize', psize);


%% panel (d) inset chi_het correction

insetD = axes('Position', [left+width+0.05, bottom+0.2, width*0.48, height*0.47]);

for rp=1:length(dims)-1
    dim = dims{rp};
    ndim = str2double(dim(1:end-1));
    fts = fitdata(:,1) == ndim;
    phihats = fitdata(fts, 2)/ndim;
    csize = exp(fitdata(fts, 6));
    typcsize = exp(fitdata(fts, 6) + fitdata(fts, 7).^2/2);
    typhet = (exp(fitdata(fts, 7).^2) - 1 ) .* exp(2*fitdata(fts, 6) + fitdata(fts, 7).^2);
    ybase = interp1(mftchhet(:,1), mftchhet(:,2), csize*ndim);
    ys = (ndim*typhet./ (typcsize.*typcsize) - ybase)./ybase;
    plot(csize*ndim, ndim*ys, 'x:', 'markersize', msize-1, 'color', colors(rp,:) );
    if 1==rp, hold on, end
end

xlim([1e-3, 1e0]);
ylim([6e-1, 2e1]);
xticks([1e-3, 1e-2, 1e-1, 1e0]);
yticks([  1,4, 10]);
set(gca, 'xscale', 'log', 'yscale', 'log');
%xlabel('$\hat\Delta$', 'interpreter', 'latex');
ylabel("$d \times \delta \hat\chi'^{\mathrm{S}}_\mathrm{het}$", 'interpreter', 'latex');
set(gca, 'YAxisLocation', 'right');

set(gca, 'fontname', 'times new roman', 'fontsize', psize-2);



%% labels

annotation('textbox',[0.075, 0.792, 0.1, 0.1], 'String', '(a)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.49, 0.792, 0.1, 0.1], 'String', '(b)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

annotation('textbox',[0.075, 0.405, 0.1, 0.1], 'String', '(c)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);
annotation('textbox',[0.49, 0.405, 0.1, 0.1], 'String', '(d)', 'EdgeColor','none', 'Fontname', 'Times New Roman','fontsize',psize);

print('../figures/rlgmf_chi.eps', '-depsc');