expdir    = '/scratch/slim/zfang/Result/UQWRI/BGModel/Paper/Test1_SmoothOnSlow';
figdir    = '/home/slim/zfang/Project/Paper/Journal/SLIM.Papers.UQofWRI/fig/BG';
figdir2   = '/scratch/slim/zfang/FigureResult/SEG/2016/UQWRI/BG';
vfile     = '/scratch/slim/zfang/Result/UQWRI/BGModel/Paper/Test1_SmoothOnSlow/m_10.mat';
vtfile    = '/scratch/slim/zfang/Model/BG2D/BG450kmTrue.mat';
vifile    = '/scratch/slim/zfang/Model/BG2D/Ini_Vel.mat';
datafile1  = '/scratch/slim/zfang/Data/BG2D/DataForUQPaperBG/BGDataNoise_31.mat';
datafile2  = '/scratch/slim/zfang/Data/BG2D/DataForUQPaperBG/Test_Cent15.mat';
addpath('/home/slim/zfang/Project/FWIwithUncertainty/Code/Version2/mbin');
StatisticFile  = [expdir '/Result_WO_Sample.mat'];
curdir    = pwd;

position       = [200 2500; 1000 2500; 1800 2500];
xpos           = [ 1000 2500 4000];
plotdis        = 0; % flag to decide whether plot distribution
plotCI         = 0; % flag to decide whether plot confidence interval
plotmodel      = 0;
plotdata         = 1;
alpha          = 0.9;
stdp           = 0.7;
vrange         = [1.5, 6];
CIxlim         = [1 5];
frq_idx       = 3;
src_idx       = 46;
xrec            = [0:10:4500];

xsize          = 10;
ysize          = 5;
vp             = ReadAllData(vifile);
vp             = vp(:);
CIPaperSize    = [6,10.5];
DataPaperSize = [8 6];
DistPaperSize  = [8 8];

if ~exist(figdir,'dir')
    mkdir(figdir);
end

if ~exist(figdir2,'dir')
    mkdir(figdir2);
end

Cfile   = mfilename('fullpath');
Cfile   = [Cfile '.m'];
copyfile(Cfile,figdir);

Label          = 'BG';

cd(figdir);

opt.save                = 1;
opt.fdir                = figdir;
opt.fformat             = 'eps';
opt.Xlabel              = ['Lateral [m]'];
opt.Zlabel              = ['Depth [m]'];
opt.textfontsize        = 20;
opt.colorbartitle       = 'km/s';
opt.PaperSize           = [10 4];
xsize = 8;
ysize = 4;
Fcaxis = [1.5,4.5];
Fcaxisstd = [0.01 0.5];


load(StatisticFile);

if plotmodel > 0
% Plot std
STDv                    = reshape(STDv,n);
opt.fname               = [Label 'STD'];
opt.Title               = 'Standard Deviation';
opt.fcaxis              = Fcaxisstd;
optstd                  = opt;
optstd.fcolormap        = redblue;
plot_2D_model(STDv,n,d,o,optstd);
scommand                = ['unset DYLD_LIBRARY_PATH;convert ' opt.fname '.eps ' opt.fname '.png']; % convert file from eps to png
[status,result]=system(scommand);


% Plot mean
opt.fname               = [Label 'MEAN'];
opt.Title               = 'Mean';
opt.fcaxis              = Fcaxis;
plot_2D_model(MEANv,n,d,o,opt);
scommand                = ['unset DYLD_LIBRARY_PATH;convert ' opt.fname '.eps ' opt.fname '.png']; % convert file from eps to png
[status,result]=system(scommand);


% Plot MAP
opt.fname               = [Label 'MAP'];
opt.Title               = 'MAP';
opt.fcaxis              = Fcaxis;
plot_2D_model(MAPv,n,d,o,opt);
scommand                = ['unset DYLD_LIBRARY_PATH;convert ' opt.fname '.eps ' opt.fname '.png']; % convert file from eps to png
[status,result]=system(scommand);

% Plot True velocity
opt.fname               = [Label 'True'];
opt.Title               = 'True model';
opt.fcaxis              = Fcaxis;
plot_2D_model(ReadAllData(vtfile),n,d,o,opt);
scommand                = ['unset DYLD_LIBRARY_PATH;convert ' opt.fname '.eps ' opt.fname '.png']; % convert file from eps to png
[status,result]=system(scommand);

% Plot Initial velocity
opt.fname               = [Label 'Initial'];
opt.Title               = 'Initial model';
opt.fcaxis              = Fcaxis;
plot_2D_model(ReadAllData(vifile),n,d,o,opt);
scommand                = ['unset DYLD_LIBRARY_PATH;convert ' opt.fname '.eps ' opt.fname '.png']; % convert file from eps to png
[status,result]=system(scommand);

end


modelsize = prod(n);

%load(filename);
MAPv  = ReadAllData(vfile);
MEANv = MAPv;

Sample = [];

for i = 1:10
       str = [expdir '/Smp_' num2str(i) '.mat'];
       load(str);
       Sample = [Sample 1./sqrt(Smp)];
end


if plotdis > 0
for i = 1:size(position, 1)
	
	% plot distribution
	k 		= sub2ind(n,(position(i,1)-o(1))/d(1)+1,(position(i,2)-o(2))/d(2)+2);
	figure; fidx    = gcf;
    set(gca, 'Fontsize', 20);
	%hist(Sample(k,:),50);
    %[YY, XX] = hist(Sample(k,:),50);
    %YY       = YY / (sum(YY) * (XX(2)-XX(1)));
    XX       = [vrange(1):0.01:vrange(2)];
    YY       = normpdf(XX,MAPv(k),STDv(k))*0.01;
    YY       = YY / (sum(YY)*0.01);
    plot(XX,YY,'linewidth',2);
    XX2      = [vrange(1):0.01:vrange(2)];
    YY2      = normpdf(XX2,vp(k),stdp)*0.01;
    YY2      = YY2 / (sum(YY2)*0.01);
    hold on;
    plot(XX2,YY2,'r','linewidth',2)
    legend('Posterior distribution','Prior distribution')
	%set(gca,'YTickLabel',sprintf('%01.2d|',get(gca,'ytick')/size(Sample,2)));
        %aaa=axes;set(aaa,'yticklabel',[]);
	xlabel('Velocity [km/s]')
	%ylabel('Distribution density')
	strtitle = ['Distribution at x = ' num2str(position(i,2)) 'm, z = ' num2str(position(i,1)) 'm'];
	strfile  = ['Distr' Label 'x' num2str(position(i,2)) 'z' num2str(position(i,1))];
	title(strtitle);
	PrintFigure(fidx,DistPaperSize,strfile);
	%keyboard
	% plot MEAN trend
	

	figure; fidx    = gcf;
    set(gca, 'Fontsize', 20);
	plot(cumsum(Sample(k,:))./[1:size(Sample,2)],'linewidth',2);
	ylabel('Velocity [km/s]')
	xlabel('Number of samples')
	strtitle = ['Trend of mean at x = ' num2str(position(i,2)) 'm, z = ' num2str(position(i,1)) 'm'];
	strfile  = [figdir '/MEANtrend' Label 'x' num2str(position(i,2)) 'z' num2str(position(i,1))];
	title(strtitle);
	PrintFigure(fidx,[xsize ysize],strfile);

    % plot STD trend
	
	figure; fidx    = gcf;
    set(gca, 'Fontsize', 20);
	
	jump    = floor(size(Sample,2)/100);
	for j = 1:100
		STDtrend(j) = std(Sample(k,1:jump*j));
	end
	
	plot([1:100]*jump,STDtrend,'linewidth',2);
	ylabel('standard deviation [km/s]')
	xlabel('Number of samples')
	strtitle = ['Trend of STD at x = ' num2str(position(i,2)) 'm, z = ' num2str(position(i,1)) 'm'];
	strfile  = [figdir '/STDtrend' Label 'x' num2str(position(i,2)) 'z' num2str(position(i,1))];
	title(strtitle);
	PrintFigure(fidx,[xsize ysize],strfile);	
	
end
	
end

% plot confidence interval
if plotCI > 0

vf = MEANv;
vt = ReadAllData(vtfile);
vt = vt;

I   = floor((xpos-o(2))/d(2)) + 1;

clear CIv;

for i = 1:length(I)
	for j = 1:n(1)
		k = sub2ind(n,j,I(i));
        CIv(j,i,:) = confidence_interval(Sample(k,:),alpha);
	end
end
[z x] = odn2grid(o,d,n);

for i = 1:length(I)
    figure;
    fidx = gcf;
    %h1 = plot_CI(z,vf(:,I(i)),CIv(:,i,1),CIv(:,i,2));
    CIvu = vec(CIv(:,i,1));
    CIvd = vec(CIv(:,i,2));
    fill([CIvu;flipud(CIvd)],[vec(z);flipud(vec(z))],[7 7 7]/8);		
    hold on;
    h2 = plot(vt(:,I(i)),z,'green');
    h3 = plot(vf(:,I(i)),z,'k');
    hline = findobj(gcf, 'type', 'line');
    set(hline,'LineWidth',3);
    xlabel('Velocity [km/s]','fontsize',20);
    ylabel('Depth [m]','fontsize',20);
    ylim([z(1),z(end)]);
	xlim(CIxlim);
    titlestr = ['x = ' num2str(x(I(i))) 'm'];
    title(titlestr,'fontsize',20);
    h  = legend([h2;h3],'V True','V Final');
    set(h,'fontsize',14);
    set(gca,'fontsize',20);
    set (gca,'Ydir','reverse')
    fname = [figdir '/CI_' Label 'x' num2str(x(I(i)))];
	PrintFigure(fidx, CIPaperSize, fname)
    
end

end

% for i = 1:length(I)
%
%     figure;
%     fidx = gcf;
%     h1 = plot_CI(z,vf(:,I(i)),CIv(:,i,1),CIv(:,i,2));
%     hold on;
%     h2 = plot(vt(:,I(i)),z,'green');
%     h3 = plot(vf(:,I(i)),z,'k');
%     hline = findobj(gcf, 'type', 'line');
%     set(hline,'LineWidth',3);
%     xlabel('Velocity [km/s]','fontsize',20);
%     ylabel('Depth [m]','fontsize',20);
%     ylim([z(1),z(end)]);
%     xlim(CIxlim);
%     titlestr = ['x = ' num2str(x(I(i))) 'm'];
%     title(titlestr,'fontsize',20);
%     h  = legend([h1(1);h2;h3],'Confidence Interval','V True','V Final');
%     set(h,'fontsize',14);
%     set(gca,'fontsize',20);
%     fname = [figdir '/CI_' Label 'x' num2str(x(I(i)))];
%     PrintFigure(fidx, CIPaperSize, fname)
%
% end

if plotdata > 0
    Dt = ReadAllData(datafile2);
    Dn = ReadAllData(datafile1);
    Dt  = squeeze(Dt);
    Dn = squeeze(Dn);
    dt  = Dt(:,src_idx,frq_idx);
    dn = Dn(:,src_idx,frq_idx);
    dt = dt(:);dn = dn(:);
    figure;
    fidx = gcf;plot(xrec, real([dn, dt]),'linewidth', 2); xlabel('Receiver [m]', 'fontsize', 20); ylabel('Real part','fontsize',20);
    h = legend('Data w/ noise', 'Data w/o noise'); set(h,'fontsize',14);
    set(gca,'fontsize',20);xlim([0 4500])
    fname = [figdir '/Data_' Label 'src' num2str(src_idx) 'frq' num2str(frq_idx)];
    PrintFigure(fidx,DataPaperSize,fname);

    
end

cd(curdir)

copyfile(figdir,figdir2);







