d0=fileparts([pwd,filesep]);
d0=d0(1:strfind(d0,'trunk')+4);
force = true;
root = 'DOnumPap';

files = dir(sprintf('%s/Save/%s/LockStoch_*',d0,root));

for iii = 1:length(files)
    load(sprintf('%s/Save/%s/%s/param.mat',d0,root,files(iii).name)) ;
    app = rmfield(app,'savemat');
    file2 = dir(sprintf('%s/Save/%s/%s/0*.mat',d0,root,files(iii).name));
    k = 0;
    LockExchangePlots_num
    for jjj = length(file2):-floor(length(file2)/2):1
        %%
        if ~exist(sprintf('%s/Save/%s/%s/%05d.png',d0,root,...
                files(iii).name,jjj),'file') || force
            clf , close all
            load(sprintf('%s/Save/%s/%s/%s',d0,root,files(iii).name,file2(jjj).name)) ;
            LockExchangePlots_num
            print('-dpng','-r400',sprintf('%s/Save/%s/%s/%05d.png',d0,root,files(iii).name,jjj));
%             print('-depsc2','-painters',sprintf('%s/Save/%s/%s/%05d.eps',d0,root,files(iii).name,jjj));
            fprintf('%s/Save/%s/%s/%05d.png\n',d0,root,files(iii).name,jjj)
            k=1;
        end
    end
end


%% Evolution of Error Plots
file2 = dir(sprintf('%s/Save/%s/%s/0*.mat',d0,root,files(1).name));
Error = zeros(length(file2), length(files));
Error_Noise = zeros(length(file2), length(files));
Error_Mean = zeros(length(file2), length(files));
for iii = 1:length(files)
    load(sprintf('%s/Save/%s/%s/param.mat',d0,root,files(iii).name)) ;
    app = rmfield(app,'savemat');
    file2 = dir(sprintf('%s/Save/%s/%s/0*.mat',d0,root,files(iii).name));
    k = 0;
    if ~exist(sprintf('%s/Save/%s/%s/evolve.png',d0,root,...
            files(iii).name),'file') || force
        clf , close all
        for jjj = 1:length(file2)
            meanu = 0;
            meanv=0;
            meanrho=0;
            for i = 0:length(app.dref)-1
                %%
                %Load the DO run
                load(sprintf('%s/Save/%s/%s/%s',d0,root,files(iii).name,file2(jjj).name)) ;
                
                %Load the Deterministic run
                det = load(sprintf('%s/Save/%s%d/%05d.mat',d0,app.ref,app.dref(i + 1),...
                    k/app.PlotIntrvl),'rho','u','v');
                meanu = meanu + det.u/4;
                meanv = meanv + det.v/4;
                meanrho = meanrho + det.rho/4;
% app.ref_noise = app.ref; %Temporary Hack
                %Load the Noisy Deterministic run
% % %                 det_noise = load(sprintf('%s/Save/%s%d/%05d.mat',d0,app.ref_noise,app.dref(i + 1),...
% % %                     k/app.PlotIntrvl/2),'rho','u','v');
                
                %Realization
                %Find coefficients for actual projections
                YYT(i+1,:) =  (det.rho - rho)' * rhoi * dx * dy; %(as a check)
                
                tmp = rhoi * YYt(i+1,:)' + rho - det.rho;
                tmp1 = ui*YYt(i+1,:)'+u - det.u;
                tmp2 = vi*YYt(i+1,:)'+v - det.v;
                nrm = sqrt((det.rho'*det.rho + det.u'*det.u + det.v'*det.v)*dx*dy);
                Error(jjj,iii) =Error(jjj,iii) + sqrt((tmp' * tmp + tmp1' * tmp1 + tmp2' * tmp2)*dx*dy) / nrm;
                
% % %                 tmp = det_noise.rho - det.rho;
% % %                 tmp1 = det_noise.u - det.u;
% % %                 tmp2 = det_noise.v - det.v;
% % %                 nrm = sqrt((det.rho'*det.rho + det.u'*det.u + det.v'*det.v)*dx*dy);
% % %                 Error_Noise(jjj,iii) = Error_Noise(jjj,iii) + sqrt((tmp' * tmp + tmp1' * tmp1 + tmp2' * tmp2)*dx*dy) / nrm;
                
                YYtTime(1:4, 1:3, jjj) = YYt(1:4,1:3);
            end
            tmp = rho - meanrho;
            tmp1 = u - meanu;
            tmp2 = v - meanv;
            Error_Mean (jjj, iii) = sqrt((tmp' * tmp + tmp1' * tmp1 + tmp2' * tmp2)*dx*dy)...
                / sqrt((meanrho' * meanrho + meanu' * meanu + meanv' * meanv)*dx*dy); 
        end

        plot(linspace(0,T,length(Error)),Error(:,iii)/4,'linewidth',1)
        hold on
        plot(linspace(0,T,length(Error)),Error_Mean(:,iii),'linewidth',2)
        
        xlabel('Time'), ylabel('||E||_2 / ||sol||_2')
        print('-dpng','-r400',sprintf('%s/Save/%s/%s/evolve.png',d0,root,files(iii).name));
        fprintf('%s/Save/%s/%s/evolve.png\n',d0,root,files(iii).name)
        k=1;
        clf, close all
        symbols = {'b-','g--','r-.','m-'};
        LW = [2,1,1,1];
        for i = 1:size(YYt,2)
            subplot(4,4,i)
            for j = 1:length(app.dref)
                %We need to fix the signs of the coefficients so that they
                %do not jump around
                k=size(YYtTime,3);
                tmpcheck = (YYtTime(j,i,k));
                for k=size(YYtTime,3):-1:1
                    YYtTime(j,i,k) = YYtTime(j,i,k) * sign(YYtTime(j,i,k) * tmpcheck);
                    tmpcheck = (YYtTime(j,i,k));
                end
                plot(linspace(0,T,length(Error)),squeeze(YYtTime(j,i,:)),symbols{j},'LineWidth',LW(j))
                hold on
            end
            xlabel('Time'), ylabel('Y_i')
            axis square
        end
        print('-dpng','-r400',sprintf('%s/Save/%s/%s/YYt_evolve.png',d0,root,files(iii).name));
        print('-depsc2','-painters',sprintf('%s/Save/%s/%s/YYt_evolve.eps',d0,root,files(iii).name));

        fprintf('%s/Save/%s/%s/YYt_evolve.png\n',d0,root,files(iii).name)

    end
end
%%
clf, close all
subplot(211)
symbols = {'r-.','g--','b-'};
for i = 2:3
    plot(linspace(0,T,length(Error)),Error(:,i)/4,symbols{i},'LineWidth',1)
    hold on
    plot(linspace(0,T,length(Error)),Error_Mean(:,i),symbols{i},'LineWidth',2)    
end
% % % symbols = {'.g-.','.r--','.b-'};
% % % for i = 2:length(files)
% % %     plot(linspace(0,T,length(Error)),Error_Noise(:,i),symbols{i},'LineWidth',1)
% % %     hold on
% % % end

xlabel('Time'), ylabel('||E||_2 / ||sol||_2')
% legend('10^7','1.25e6','4e5','Location','NorthWest')
legend('10^7','1.25e6','4e4','Location','NorthWest')
print('-dpng','-r300',sprintf('%s/Save/%s/%s/../LE_TOTevolve.png',d0,root,files(iii).name));
print('-depsc2',sprintf('%s/Save/%s/%s/../LE_TOTevolve.eps',d0,root,files(iii).name));
fprintf('%s/Save/%s/%s/TOTevolve.png\n',d0,root,files(iii).name)

%%
subplot(4,4,1)
plot(repmat([app.dref/100-0.8]'/2,1,2)',[0,0,0,0;1,1,1,1],'b','linewidth',1)
xlim([-0.25,0.25]/2)
axis square
print('-dpng','-r400',sprintf('%s/Save/%s/%s/../LEPDF_init.png',d0,root,files(iii).name));