% MSEplot
% go
% go2 with Sample Size of 250!
%
%%% msePLOT
fwd = 30;
A1  = zeros(size(anlyc,1),size(anlyc,2),fwd+1);
AC  = zeros(size(anlyc,1),size(anlyc,2),fwd+1);
A1(:,:,1) = anly1(:,:,2);
AC(:,:,1) = anlyc(:,:,2);
RS  = TRU(:,11);
warning off;
RS  = forward_1d(RS,fwd,1);


tic; k = 1;
for i=1:size(anly1,2);
    A1(:,i,:) = squeeze(forward_1d(A1(:,i,1),fwd,0));
    AC(:,i,:) = squeeze(forward_1d(AC(:,i,1),fwd,0));
    if(k==25); disp(toc); k=0; end
    k = k+1;
end
warning on;

RS = RS(1:100,:);
A1 = A1(1:100,:,:);
AC = AC(1:100,:,:);
xx = msh;
nc = size(A1,2)

f1 = quantile(fore1(1:100,:,1),[.95,.05],2);
a1 = quantile(A1(:,:,1),[.95,.05],2);
ac = quantile(AC(:,:,1),[.95,.05],2);

plot(msh,f1); ax = axis;
plot(msh,f1(:,1),'b',msh,a1(:,1),'r:',msh,ac(:,1),'g--')
axis(ax)

hold on; plot(1:100,A1(:,:,1),'y:'); hold off
hold on; plot(1:100,AC(:,:,1),'r:'); hold off
hold on; plot(1:100,RS(:,1),'k'); hold off


biasC = zeros(size(AC,1),length(fwd)); 
bias1 = biasC;
mseC  = biasC;  mse1  = bias1;
madC  = biasC;  mad1  = bias1;
for i = 1:fwd
    dum = repmat(RS(:,i),[1,Nens]);
    biasC(:,i) = mean( AC(:,:,i)-dum,2);
    mseC(:,i)  = sqrt(mean((AC(:,:,i)-dum).^2,2));
    madC(:,i)  = median(abs(AC(:,:,i)-dum),2);
    bias1(:,i) = mean( A1(:,:,i)-dum,2);
    mse1(:,i)  = sqrt(mean((A1(:,:,i)-dum).^2,2));
    mad1(:,i)  = median(abs(A1(:,:,i)-dum),2);
end

j = 5;
plot(msh,mse1(:,j),'b',msh,mseC(:,j),'r')
plot(msh,bias1(:,j),'b',msh,biasC(:,j),'r')

plot(msh, sqrt(mse1(:,5)),'b')

for j = 1:fwd
    plot(msh,A1(:,:,j),'b:',msh,AC(:,:,j),'r',msh,RS(:,j),'k')
    title(sprintf('%d',j))
    pause(.02)
end

imagesc(mse1(:,1:10)'); nrm; colorbar
imagesc(mseC(:,1:10)'); nrm; colorbar

r = (1:fwd) +10;
[xx,yy] = meshgrid(msh,r);
%plot3(xx,yy,squeeze(mse1(:,r)./mseC(:,r)))

figure
subplot(3,1,1)
  plot(msh,fore1(1:100,:,1),'g')
  line(msh,TR(:,1),'Color','red','LineWidth',2)
  ylabel('Temperature');  axis([0,.65,-100,2400])
  title(' (a) Forecast')
subplot(3,1,2)
  plot(msh,anly1(1:100,:,2),'g')
  line(msh,TR(:,1),'Color','red','LineWidth',2)
  hold on; plot(msh(spc(1:10)),dat(1:10,1),'k^'); hold off
  ylabel('Temperature')
  title(' (b) Analysis'); axis([0,.65,-100,2400])
subplot(3,1,3)
  plot(msh,anlyc(1:100,:,2),'g')
  line(msh,TR(:,1),'Color','red','LineWidth',2)
  hold on; plot(msh(spc(1:10)),dat(1:10,1),'k^'); hold off
  ylabel('Temperature');  axis([0,.65,-100,2400])
  title(' (c) Penalized Analysis')
  xlabel('Spatial Location'); 
  


createMESHfigure(yy',xx',log(mse1(:,r)./mseC(:,r)))
    xlabel('Time Iteration')
    ylabel('Spatial Location')
    zlabel('log( MSE_A/MSE_C )')

createMESHfigure(yy',xx',TR(:,r))
    xlabel('Time Iteration')
    ylabel('Spatial Location')
    zlabel('Temperature')
    
createMESHfigure(yy',xx',log(abs(bias1(:,r))./abs(biasC(:,r))))
%    xlabel('Time Iteration')
%    ylabel('Spatial Location')
%    zlabel('log( bias_A/bias_C )')


j = 2;
A31 = anly1(:,:,j+1);
A3C = anlyc(:,:,j+1);
plot(msh,A31(1:100,:),'b',msh,A3C(1:100,:),'r:')
fwd31 = zeros([size(A31),fwd+1]);
fwd3c = fwd31;


warning off;
tic; k = 1;
for i=1:Nens;
    fwd31(:,i,:) = squeeze(forward_1d(anly1(:,i,j+1),fwd,0));
    fwd3c(:,i,:) = squeeze(forward_1d(anlyc(:,i,j+1),fwd,0));
    if(k==25); disp([i,round(toc)]); k=0; end
    k = k+1;
end
warning on;