function [iter,Jo,Jb,gradJ]=cost_function(fname,normalize);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2002 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [iter,Jo,Jb,GradJ]=cost_function(fname,normalize)                %
%                                                                           %
% This function reads and plots 4DVAR cost function.                        %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    fname       NetCDF initial file name (character string).               %
%    normalize   Normalize cost function by first minimization value.       %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    iter        Iteration (matrix).                                        %
%    Jo          Observation cost function (matrix).                        %
%    Jb          Background cost function (matrix).                         %
%    gradJ       total cost function gradient norm (matrix).                %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin < 2),
  normalize=0;
end,

% Inquire about variables.

[vname,nvars]=nc_vname(fname);

got_gradJ=0;
got_Jb=0;
got_Jo=0;
got_Jmin=0;

for n=1:nvars
  name=deblank(vname(n,:));
  switch name
    case 'cost_gradient'
      got_gradJ=1;
    case 'back_function'
      got_Jb=1;
    case 'cost_function'
      got_Jo=1;
    case 'Jmin'
      got_Jmin=1;
  end,
end,

%  Read in cost function data.

spval=1.0e+25;

if (got_gradJ),
  gradJ=nc_read(fname,'cost_gradient');
  [mgrd,ngrd]=size(gradJ);
end,
if (got_Jb),
  Jb=nc_read(fname,'back_function');
  [mbck,nbck]=size(Jb);
end,
if (got_Jo),
  Jo=nc_read(fname,'cost_function');
  [mobs,nobs]=size(Jo);
end,
if (got_Jmin),
  Jmin=nc_read(fname,'Jmin');
end,

%  Set labels.

if (got_Jo),
  if (length(Jo(:,1)) > 4),
    state={'Total','zeta','ubar','vbar','u','v','temp','salt'};
  else,
    state={'Total','zeta','ubar','vbar'};
  end,
end,

%  Remove fill values and reshape matrix size.

J_state=state;
J_index(1:mbck)=1;
ic=1;
for m=2:mobs,
  if (Jo(m,1) > 0 & Jo(m,1) < spval),
    ic=ic+1;
  else,
    J_index(m)=NaN;
  end,
end,
mobs=ic;

Jo=Jo(:);
i=find(Jo > spval);
if (~isempty(i)),
  Jo(i)=NaN;
end,
i=find(Jo == 0);
if (~isempty(i)),
  Jo(i)=NaN;
end,

i=find(isnan(Jo));
if (~isempty(i)),
  Jo(i)=[];
  nobs=fix(length(Jo)/mobs);
  Jo=reshape(Jo,mobs,nobs);
  if (got_Jb),
    JB=Jb(:);
    JB(i)=[];
    JB=reshape(JB,mobs,nobs);
  end,
end,

i=find(isnan(J_index));
if (~isempty(i)),
  J_state(i)=[];
end,


if (got_Jb),
  Jb_state=state;
  Jb_index(1:mbck)=1;
  Jb(:,nobs+1:end)=NaN;
  ic=1;
  for m=2:mbck,
    if (Jb(m,2) == 0),
      Jb(m,:)=NaN;
      Jb_index(m)=NaN;
    elseif (Jb(m,2) > spval),
      Jb(m,:)=NaN;
      Jb_index(m)=NaN;
    else,
      ic=ic+1;
    end,
  end,
  mbck=ic;
  nbck=nobs;
  Jb=Jb(:);
  [i]=find(isnan(Jb));
  if (~isempty(i)),
    Jb(i)=[];
  end,
  Jb=reshape(Jb,mbck,nbck);

  i=find(isnan(Jb_index));
  if (~isempty(i)),
    Jb_state(i)=[];
  end,
end,

if (got_gradJ)
  gradJ_state=state;
  gradJ_index(1:mgrd)=1;
  gradJ(:,nobs+1:end)=NaN;
  ic=1;
  for m=2:mgrd,
    if (gradJ(m,2) == 0),
      gradJ(m,:)=NaN;
      gradJ_index(m)=NaN;
    elseif (gradJ(m,2) > spval),
      gradJ(m,:)=NaN;
      gradJ_index(m)=NaN;
    else,
      ic=ic+1;
    end,
  end,
  mgrd=ic;
  ngrd=nobs;
  gradJ=gradJ(:);
  [i]=find(isnan(gradJ));
  if (~isempty(i)),
    gradJ(i)=[];
  end,
  gradJ=reshape(gradJ,mbck,nbck);

  i=find(isnan(gradJ_index));
  if (~isempty(i)),
    gradJ_state(i)=[];
  end,
end,

%  Set default color order.

Corder=[0.00 0.00 0.00; ...  % black
        1.00 0.00 0.00; ...  % red
        0.00 1.00 1.00; ...  % cyan
        0.00 1.00 0.00; ...  % green
        0.00 0.00 1.00; ...  % blue
        1.00 0.00 1.00; ...  % magenta
        1.00 0.75 0.00; ...  % orange
        1.00 0.50 0.25];     % pink

set(0,'DefaultAxesColorOrder',Corder);

if (got_Jb),
  J=JB+Jo;
else,
  J=Jo;
end,
J=J';

% Plot observation cost function.

if (got_Jo),
  figure;

  x=1:nobs;
  iter=repmat(x,[mobs,1]);
  iter=iter';

  Jo=Jo';
  if (normalize),
    for i=1:mobs,
      Jo(:,i)=Jo(:,i)./J(2,i);
    end,
  end,

  set(gcf,'Units','Normalized',...
          'Position',[0.2 0.1 0.6 0.8],...
	  'PaperUnits','Normalized',...
	  'PaperPosition',[0.2 0.1 0.6 0.8]);

  semilogy(iter,Jo,'.-')
  legend(J_state,1);
  xlabel('Iteration');
  title('4DVAR Observation Cost Function, J_o');
  grid on;
  print -dpsc cost.ps
end,

%  Plot background cost function

if (got_Jb),
  figure;

  x=1:nbck;
  iter=repmat(x,[mbck,1]);
  iter=iter';

  Jb=Jb';
  if (normalize & (mbck == mobs)),
    for i=1:mbck,
      Jb(:,i)=Jb(:,i)./J(2,i);
    end,
  end,

  set(gcf,'Units','Normalized',...
          'Position',[0.2 0.1 0.6 0.8],...
	  'PaperUnits','Normalized',...
	  'PaperPosition',[0.2 0.1 0.6 0.8]);

  semilogy(iter,Jb,'.-')
  legend(Jb_state,'Location','SouthEast');
  xlabel('Iteration');
  title('4DVAR Background Cost Function, J_b');
  grid on;
  print -dpsc -append cost.ps
end,

% Plot cost function gradient norm.

if (got_gradJ),
  figure;

  x=1:ngrd;
  iter=repmat(x,[mgrd,1]);
  iter=iter';

  gradJ=gradJ';
  if (normalize),
    for i=1:mgrd,
      gradJ(:,i)=gradJ(:,i)./gradJ(1,i);
    end,
  end,

  set(gcf,'Units','Normalized',...
          'Position',[0.2 0.1 0.6 0.8],...
	  'PaperUnits','Normalized',...
	  'PaperPosition',[0.2 0.1 0.6 0.8]);

  semilogy(iter,gradJ,'.-')
  legend(gradJ_state,1);
  xlabel('Iteration');
  title('4DVAR Cost Function Gradient Norm');
  grid on;
  print -dpsc -append cost.ps
end,

% Plot optimality property.

if (got_Jmin),
  figure;

  x=1:length(Jmin);
  y=ones(length(Jmin));

  set(gcf,'Units','Normalized',...
          'Position',[0.2 0.1 0.6 0.8],...
	  'PaperUnits','Normalized',...
	  'PaperPosition',[0.2 0.1 0.6 0.8]);

  plot(x,Jmin,'.-',x,y,'r-')
  legend('Estimated','Theoretical');
  xlabel('Iteration');
  title('Optimal, Normalized Cost Function Minimun');
  grid on;
  print -dpsc -append cost.ps
end,


return

