function varargout=mleos(Hx,Gx,thini,params,blurs)
% [thhat,logli,thini,scl,params,eflag,oput,grd,hes]=...
%          MLEOS(Hx,Gx,thini,params,blurs)
%
% Performs a maximum-likelihood estimation for UNCORRELATED loads as in
% Olhede & Simons by unconstrained minimization using Matlab's FMINUNC.
%
% INPUT:
%
% Hx       Real matrix with surface and subsurface topography [m m], 
% Gx       Real vector with Bouguer gravity anomaly data [mgal]
%          ... either from direct observation or from SIMULOS...
%          ... with the geographical coordinates linearly unwrapped...
% thini    A starting guess for the parameter vector with elements:
%          [D f2 s2 nu rho], in Nm, and "nothing", see SIMULOS
% params   A parameter structure with constants assumed known, see SIMULOS
%          [DEL g z2 dydx NyNx] in kg/m^3, kg/m^3, m/s^2, m, m and "nothing"
% blurs    1 Blur likelihood using the Fejer window [default]
%          0 Don't blur likelihood using the Fejer window
%
% OUTPUT:
%
% thhat    The maximum-likelihood estimate of the vector with elements:
%          [D f2 s2 nu rho], in Nm, and "nothing", see SIMULOS
% logli    The maximized value of the likelihood
% thini    The starting guess used in the optimization
% scl      The scaling applied as part of the optimization procedure
% params   The known constants used inside, see above
% eflag    The exit flag of the optimization procedure [bad if 0]
% oput     The output structure of the optimization procedure
% grd      The gradient at the estimate
% hes      The Hessian at the estimate
%
% EXAMPLE:
%
%% Perform a series of simulations
% mleos('demo1',100)
%% Statistical study of a series of simulations
% mleos('demo2','10-May-2011')
%% Admittance/coherence study of a series of simulations
% mleos('demo3','10-May-2011')
%% Covariance study of a series of simulations
% mleos('demo4','10-May-2011')
%
% Last modified by fjsimons-at-alum.mit.edu, 03/16/2012

if ~isstr(Hx)
  defval('blurs',1)

  % The necessary strings for formatting
  str0='%15s';
  str1='%12.0e';
  str2='%12.3g';
  
  % Get or supply the needed parameters
  fields={'DEL','g','z2','dydx','NyNx'};
  defstruct('params',fields,...
	    {[2670 630],9.81,35000,[20 20]*1e3,sqrt(length(Gx))*[1 1]});
  
  % Extract the variables explicitly from this structure
  for ind=1:length(fields)
    eval(sprintf('%s=params.(fields{ind});',fields{ind}))
  end
  
  % The gravitational constant (in m^3/kg/s2)
  G=fralmanac('GravCst');

  % Being extra careful or not?
  defval('xver',1)
  
  % The random guesses for the parameters, see also SIMULOS
  aguess=[1e24 0.8 0.0025 2 2e4];
  % Variations of this will need to be incorporated into this routine
  aguess=[7e22 0.4 0.0025 2 2e4];
  
  % The number of parameters to solve for
  np=length(aguess);
  % The number of unique entries in an np*np symmetric matrix
  npp=np*(np+1)/2;
  
  % Scale the parameters by this factor; stick with it throughout
  scl=10.^round(log10(aguess));
  aguess=aguess./scl;
  
  % Perturb the guesses but keep them all positive if possible
  nperturb=0.25;
  % So not all the initialization points are the same!!
  defval('thini',max((1+nperturb*[randn(size(aguess))+3]).*aguess,aguess))
  
  % But if you supplied your own guess it needs to be scaled!
  
  % Talk to me
  disp(sprintf(sprintf('%s : %s ',str0,repmat(str1,size(scl))),...
	       'Scaling',scl))
  disp(sprintf(sprintf('%s : %s ',str0,repmat(str2,size(thini))),...
	       'Starting theta',thini.*scl))
  
  % Create the appropriate wavenumber axis
  k=knum2(NyNx,[(NyNx(1)-1)*dydx(1) (NyNx(2)-1)*dydx(2)]);

  defval('taper',0)
  if taper==1
    % Were going to want to make a 2D taper - any taper
    disp('MLEOS with TAPERING, DO NOT DO THIS YET')
    NW=2;
    E1=dpss(NyNx(1),NW,1);
    E2=dpss(NyNx(2),NW,1);
    Tx=repmat(E1,1,NyNx(2)).*repmat(E2',NyNx(1),1);
    % But should still watch out for the spectral gain I suppose, this isn't
    % done just yet, although it appears already properly normalized
    % However, this looks better now, doesn't it?
    Tx=Tx*sqrt(prod(NyNx));
    % Not doing anything still amounts to saying Tx=1
  else
    Tx=1;
  end
    
  % Turn the observation vector to the spectral domain
  Hk(:,1)=tospec(Tx(:).*Hx(:,1),NyNx);
  Gk     =tospec(Tx(:).*Gx     ,NyNx);
  Hk(:,2)=Gk.*exp(k(:).*z2)/2/pi/G/DEL(2);

  if xver==1
    % This is reliant on the gravity data to be amenable to deconvolution,
    % while we could fake it in the simulations by working with Hx. Check
    % quickly, and note that there are roundoff errors right away!
    % Is the normalization right? I recently absorbed this into TOSPEC.
    difer([tospec(Tx(:).*Hx(:,2),NyNx)-Hk(:,2)]/length(Hk),8,[],NaN)
    % Should also compare this with what actually can come out of SIMULOS
    % itself, although with real data of course we don't have this.
  end
  
  NN=200;
  % And now get going with the likelihood using Hk(:,1:2) or [Hk(:,1) Gk]
  % [ off | iter | iter-detailed | notify | notify-detailed | final |
  % final-detailed ] 
  options=optimset('GradObj','off','Display','iter-detailed',...
		   'TolFun',1e-11,'TolX',1e-11,'MaxIter',NN,...
		   		   'LargeScale','off');
  % The 'LargeScale' option goes straight for the line search when the
  % gradientis NOT being supplied.

  % Set the parallel option to (never) use it for the actual optimization
  % Doesn't seem to do much when we supply our own gradient
  options.UseParallel='always';

  if xver==1 && blurs~=1
    options.DerivativeCheck='on';
  end

  % And find the MLE! Work on scaled parameters
  try
    [thhat,logli,eflag,oput,grd,hes]=...
	fminunc(@(theta) loglios(theta,params,Hk,k,scl,blurs),...
		thini,options);
  catch
    varargout=cellnan(nargout,1,1);
    return
  end
  
  % When you're done cold compared grd and hes with our own ?
    
  % Now if this was a single estimate that we are making we'd use COVTHOS
  % here to compute its covariance matrix. In 'demo1' we make sure that
  % this makes sense.
  
  % Talk!
  disp(sprintf(sprintf('%s : %s ',str0,repmat(str2,size(thhat))),...
	       'Estimated theta',thhat.*scl))

  % Generate output as needed
  varns={thhat,logli,thini,scl,params,eflag,oput,grd,hes};
  varargout=varns(1:nargout);
elseif strcmp(Hx,'demo1')
  %if ~matlabpool('size')
    %matlabpool open
  %end
  % If you run this again on the same date, we'll just add to THINI and
  % THHAT but you should start WITH A BLANK THZERO. See 'demo2'
  % How many simulations: the second argument after the demo id
  defval('Gx',[]);
  N=Gx;
  defval('N',1000)
  more off

  % The number of parameters to solve for
  np=6;
  % The number of unique entries in an np*np symmetric matrix
  npp=np*(np+1)/2;

  % Ouput files, in parallel might be a jumble
  % The thruth and the theoretical covariances
  fid0=fopen(sprintf('thzro_%s',date),'w');
  % The estimates
  fid1=fopen(sprintf('thhat_%s',date),'a+');
  % The initial guesses
  fid2=fopen(sprintf('thini_%s',date),'a+');
  % The collected optimization diagnostics, replicated some above
  fid3=fopen(sprintf('diagn_%s',date),'a+');

  % Output formating for the estimation parameters
  fmt1=['%12.6e ' repmat('%12.6f ',1,np-1) '\n'];
  % Output formating for the simulation parameters
  fmt2='%i %i %f %i %f %f %i %i\n';
  % For the time, exit flag, iterations 
  fmta='%3i %3i %3i\n';
  % For the likelihood and first-order optimality
  fmtb='%15.8e %15.8e\n';
  % For the scale
  fmtc=[repmat('%15.0e ',1,np) '\n'];
  % For the score
  fmtd=[repmat('%15.8e ',1,np) '\n']; 
  % For the Hessian, with 15 unique elements:
  fmte=repmat([repmat('%15.12f ',1,npp/3) '\n'],1,3); 
  fmt3=[fmta fmtb fmtc fmtd fmte];
  % For the unscaled covariance matrix
  fmtg=repmat([repmat('%19.12e ',1,npp/3) '\n'],1,3);

  % Do it!
  good=0; 
  % Initialize the average Hessian
  avhs=sparse(np,np);
  
  for index=1:N
    % Simulate data from the same lithosphere
    [Hx,Gx,th0,p,k]=simulos;
    % Check the dimensions of space and spectrum are right
    difer(length(Hx)-length(k(:)),[],[],NaN)
    % Initilize the THZRO file
    if index==1
      % Commit the truth to file
      fprintf(fid0,fmt1,th0);
      % Commit the parameters of the experiment to file
      fprintf(fid0,fmt2,struct2array(p));
      % Close and reopen?
    end
    % Form the maximum-likelihood estimate
    t0=clock;
    [thhat,logli,thini,scl,p,e,o,gr,hs]=mleos(Hx,Gx);
    ts=etime(clock,t0);

    % If a model was found, keep the results
    % Ignore the fact that it may be at the maximum number of iterations
    % e=1
    itmin=10;
    optmin=1;
    %if e>0 && isreal([logli gr' full(hs(:))']) && all(thhat>0) && all(~isnan(thhat))
    if isreal([logli gr' full(hs(:))']) ...
	  && all(thhat>0) ...
	  && all(~isnan(thhat)) ...
	  && o.iterations > itmin ...
	  && o.firstorderopt < optmin
      good=good+1;
      % Build the average of the Hessians for printout later
      avhs=avhs+hs;
      % Reapply the scaling before writing it out
      fprintf(fid1,fmt1,thhat.*scl);
      fprintf(fid2,fmt1,thini.*scl);

      % Print the optimization diagnostics to a different file
      hessi=tril(full(hs));
      
      % First line: The estimate, scaled back to proper units
      fprintf(fid3,fmt1,thhat.*scl);
      % First line: The initial guess, scaled back to proper units
      fprintf(fid3,fmt1,thini.*scl);
      % Third line: The diagnostics, likelihood, scales, and derivatives
      fprintf(fid3,fmt3,...
	      round(ts),e,o.iterations,...
	      logli,o.firstorderopt,...
	      scl,gr,hessi(~~hessi));
    end
  end

  % If there was any success at all, finalize the THZERO file 
  % If for some reason this didn't end well, simply rerun a couple of
  % iterations and steal the last bits. The average Hessian is likely
  % junk anyway, and you just need the covariance and Fisher which don't
  % change. See the note below also.
  
  % Initialize if all you want is to close the file
  defval('scl',nan(np,1))
  defval('th0',nan(np,1))
  if N==0; [Hx,Gx,th0,p,k]=simulos; good=1; end
  
  if good>=1
    % Print the single scaling of the theoretical values that was used 
    fprintf(fid0,fmtc,scl);

    % This is the average of the Hessians, should be close to the Fisher
    avhs=full(avhs/good);
    avlin=tril(avhs); avlin=avlin(~~avlin)';
    % Print the average of the Hessians from the simulations
    fprintf(fid0,fmte,avlin);
    
    % Now compute the theoretical covariance and scaled Fisher
    % This we don't really need any other time, it's just for use to be
    % able to compare to the average of the Hessians in this simulation
    sclth0=10.^round(log10(th0));
    [covth,covlin,Fisher,Flin]=covthos(th0./sclth0,p,k,sclth0);
    fprintf(fid0,fmte,Flin);
    
    % Now print the theoretical covariance to file also
    % Of course when we don't have the truth we'll build the covariance
    % from the single estimate that we have just obtained. This
    % covariance would then be the only thing we'd have to save.
    covs=tril(covth);
    fprintf(fid0,fmtg,covs(~~covs));
    % Kind of forget what I was going to do with covlin, just let it go. 
  end
  
  % Put both of these also into the thzro file 
  fclose('all');
  
  %matlabpool close
elseif strcmp(Hx,'demo2')
  defval('Gx',[]);
  datum=Gx;
  defval('datum',date)

  % The number of parameters to solve for
  np=5;
  % The number of unique entries in an np*np symmetric matrix
  npp=np*(np+1)/2;

  % You need to be in the proper directory to do this
  f1=sprintf('thzro_%s',datum);
  
  % Then load and display what we have
  fid=fopen(f1,'r');
  % The truth
  th0=fscanf(fid,'%f',np)';
  
  % The other parameters of the experiment
  params=fscanf(fid,'%f',8)';
  % The scale used in the optimization
  sclth0=fscanf(fid,'%f',np)';
  % Don't necessarily look at THIS partial average of the Hessians
  % through the iterations as it's just of the last few iterations that
  % add cumulatively to THINI and THHAT. We are getting this again from
  % the full file diagnos. So if we have interrrupted a sequence of
  % simulations we need to provide this, even if we make it up. Or
  % rather, we would have to run one more simulation to close out this
  % file properly.
  likelyjunk=fscanf(fid,'%f',npp)';
  % But this is the scaled Fisher matrix
  Fisher=fscanf(fid,'%f',npp)';
  % And this is the theoretical covariance of the estimate
  truecov=fscanf(fid,'%f',npp)';
  % This is what we had first done wrong
  fclose(fid);
  
  % Now reassemble the covariance
  fullcov=zeros(np,np);
  fullcov(nonzeros(triu(reshape(1:np^2,np,np)')'))=truecov;
  truecov=[tril(fullcov)'+tril(fullcov)-diag(diag(fullcov))];

  % Load the initial guess 
  f2=sprintf('thini_%s',datum);
  thinis=load(f2);
  % Load the the estimates
  f3=sprintf('thhat_%s',datum);
  thhats=load(f3);
  % Load the optimization diagnostics, which should duplicate f2 and f3
  f4=sprintf('diagn_%s',datum);
  [thhat,thini,tseiter,scl,L,gam,hes]=diagnos(f4,pwd); 

  % Check that this last file captured it all and then adapt later to let
  % the individual file reads go
  difer(thhats-thhat,[],[],NaN)
  difer(thini-thinis,[],[],NaN)

  % See if we need to worry about scaling at all (if they differed,
  % although for the analysis this is totally immaterial)
  difer(sum(sum(diff(scl,1,1),1)),[],[],NaN)
  % So this observed average is to be compared with the theory
  avhs=mean(hes,1);
  
  % The necessary strings for formatting
  str0='%15s';
  str1='%12.5e ';
  str2='%12.3g';

  % Diagnostics from the simulations compared to theory
  disp(sprintf(sprintf('%s : %s ',str0,repmat(str1,size(th0))),...
	       'True theta',th0))
  disp(sprintf(sprintf('%s : %s ',str0,repmat(str1,size(th0))),...
	       'Average estimated theta',mean(thhats,1)))

  % These two should be close!
  disp(sprintf('Over %i simulations, the average Hessian and',size(thinis,1)))
  disp(sprintf(...
      'the Fisher matrix are |%4.2f|%s apart on average',...
      1/100*round(100*mean(abs([avhs-Fisher]'./Fisher'*100))),'%'))

  % Check the empirical (co)variance also
  disp(sprintf('%15s : %12.5e %12.5e %12.5e %12.5e %12.5e',...
	       'Observed standard deviation',std(thhats)))
  disp(sprintf('%15s : %i %i %i %i %i %s',...
	       'Observed standard deviation (relative)',...
	       round(100*std(thhats)./th0),'%'))
  disp(sprintf('%15s : %12.5e %12.5e %12.5e %12.5e %12.5e',...
	       'Theortcl standard deviation',sqrt(diag(truecov))))
  disp(sprintf('%15s : %i %i %i %i %i %s',...
	       'Ratio of observation to theory',...
	       round(100*std(thhats)./sqrt(diag(truecov)')),'%'))
  
  % Easiest to compare the scaled versions
  % simulcov=cov(thhats)./[sclth0(:)*sclth0(:)'];
  % scltruec=truecov./[sclth0(:)*sclth0(:)'];

  % Plot it all
  [ah,ha]=mleplos(thhats,th0,truecov);
  
  % Stick the params here somewhere so we can continue to judge
  movev(ah,-.1)
  t=supertit(ah(1:length(ah)/2),sprintf(...
      '%s %s_1 = %i ; %s_2 = %i %s ; z_2 = %i km ; %ix%i grid ; %ix%i km',...
      'MLE Simulations with',...
      '\Delta',params(1),'\Delta',params(2),'kg m^{-3}',...
      params(4)/1e3,params(7:8),params(7:8).*params(5:6)/1e3));
  movev(t,.1)
  % Young's modulus 
  defval('E',1.4e11);
  % Poisson's ratio
  defval('v',0.25);
  % Conversion to Te [m]
  Te=DtoTe(th0(1),E,v);
  
  omukm=round(mean(DtoTe(thhats(:,1),E,v)/1000));
  ostkm=round(std(DtoTe(thhats(:,1),E,v)/1000));
  disp(sprintf('observed mean %i km',omukm))
  disp(sprintf('observed stdv %3.2f km',ostkm))

  % This will come awfully close but MUST PUT IN THE MEAN!
  fakit=DtoTe(th0(1)+randn(10000,1)*sqrt(truecov(1,1)),E,v);
  stdTe=std(fakit);

  disp(sprintf('threticl std %3.2f km',stdTe/1000))

  stdf2=sqrt(truecov(2,2));
  stds2=sqrt(truecov(3,3));
  stdnu=sqrt(truecov(4,4));
  stdrho=sqrt(truecov(5,5));
  
  % We are quoting the TRUTHS and the THEORETICAL STANDARD DEVIATION with
  % which it can be known using the available data
  answ1=sprintf('E = %g',E);
  answ2=sprintf('%s = %g','\nu',v);
  answ3=sprintf('T_e = %4.1f %s %3.1f km',       round(Te/1e3*10)/10,str2mat(177),round(stdTe/1e3*10)/10);
  answ4=sprintf('f^2 = %g %s %g',                th0(2),str2mat(177),round(stdf2*1e3)/1e3);
  answ5=sprintf('%s^2 = %6.4f %s %6.4f','\sigma',th0(3),str2mat(177),stds2);
  answ6=sprintf('%s = %5.3f %s %5.3f','\nu',     th0(4),str2mat(177),stdnu);
  answ7=sprintf('%s = %i %s %i','\rho',          round(th0(5)),str2mat(177),round(stdrho));

  tt=supertit(ah(round(3*length(ah)/4)),...
	      sprintf('%s ; %s ; %s ; %s\n %s ; %s ; %s',answ1,answ2,answ3,answ4,answ5,answ6,answ7));
  movev(tt,-4)
  % Print the figure! Don't forget the degs Perl script
  figdisp([],sprintf('%s_%s',Hx,datum),[],0)
elseif strcmp(Hx,'demo3')
  defval('Gx',[]);
  datum=Gx;
  defval('datum',date)

  % The number of parameters to solve for
  np=5;
  % The number of unique entries in an np*np symmetric matrix
  npp=np*(np+1)/2;

  % You need to be in the proper directory to do this, e.g.
  % /u/fjsimons/MyPapers/OLHEDE/Simulations
  f1=sprintf('thzro_%s',datum);
  f4=sprintf('diagn_%s',datum);
  
  % Then load and display what we have
  fid=fopen(f1,'r');
  % The truth
  th0=fscanf(fid,'%f',5)';
  % The other parameters of the experiment, see SIMULOS and MLEOS
  params=fscanf(fid,'%f',8)';
  % The scale used in the optimization
  sclth0=fscanf(fid,'%f',5)';
  % Don't necessarily look at THIS partial average of the Hessians
  % through the iterations as it's just of the last few iterations that
  % add cumulatively to THINI and THHAT. We are getting this again from
  % the full file diagnos.
  likelyjunk=fscanf(fid,'%f',npp)';
  % But this is the scaled Fisher matrix
  Fisher=fscanf(fid,'%f',npp)';
  % And this is the theoretical covariance of the estimate
  truecov=fscanf(fid,'%f',npp)';
  fclose(fid);
  % Now reassemble the covariance
  fullcov=zeros(np,np);
  fullcov(nonzeros(triu(reshape(1:np^2,np,np)')'))=truecov;
  truecov=[tril(fullcov)'+tril(fullcov)-diag(diag(fullcov))];
  
  % Load the initial guess and the estimates
  f2=sprintf('thini_%s',datum);
  thinis=load(f2);
  f3=sprintf('thhat_%s',datum);
  thhats=load(f3);

  % Load the optimization diagnostics, which should duplicate f2 and f3
  [thhat,thini,tseiter,scl,L,gam,hes]=diagnos(f4,pwd); 
  
  % Check that this last file captured it all and then adapt later to let
  % the individual file reads go
  difer(thhats-thhat,[],[],NaN)
  difer(thini-thinis,[],[],NaN)

  % Plot it all: one admittance/coherence curve for every estimate
  % Well, only plot a hundred at random, how about that
  [ah,ha]=admitos(thhats(randi(length(thhats),100,1),:),th0,params);
  
  % Final cosmetics
  shrink(ah,1,2)
  fig2print(gcf,'portrait')
  t=supertit(ah,sprintf(...
      '%s %s_1 = %i ; %s_2 = %i %s \n z_2 = %i km ; %ix%i grid ; %ix%i km',...
      sprintf('%i MLE Simulations with',length(thhats(:,1))),...
      '\Delta',params(1),'\Delta',params(2),'kg m^{-3}',...
      params(4)/1e3,params(7:8),params(7:8).*params(5:6)/1e3));
  movev(t,0.4)

  % Young's modulus 
  defval('E',1.4e11);
  % Poisson's ratio
  defval('v',0.25);

  % Conversion to Te [m]
  Te=DtoTe(th0(1),E,v);
  
  % This will come awfully close but MUST PUT IN THE MEAN!
  fakit=DtoTe(th0(1)+randn(10000,1)*sqrt(truecov(1,1)),E,v);
  stdTe=std(fakit);
  stdf2=sqrt(truecov(2,2));
  stds2=sqrt(truecov(3,3));
  stdnu=sqrt(truecov(4,4));
  stdrho=sqrt(truecov(5,5));
  
  % We are quoting the TRUTHS and the THEORETICAL STANDARD DEVIATION with
  % which it can be known using the available data
  answ1=sprintf('E = %g',E);
  answ2=sprintf('%s = %g','\nu',v);
  answ3=sprintf('T_e = %4.1f %s %3.1f km',       round(Te/1e3*10)/10,str2mat(177),round(stdTe/1e3*10)/10);
  answ4=sprintf('f^2 = %g %s %g',                th0(2),str2mat(177),round(stdf2*1e3)/1e3);
  answ5=sprintf('%s^2 = %6.4f %s %6.4f','\sigma',th0(3),str2mat(177),stds2);
  answ6=sprintf('%s = %5.3f %s %5.3f','\nu',     th0(4),str2mat(177),stdnu);
  answ7=sprintf('%s = %i %s %i','\rho',          round(th0(5)),str2mat(177),round(stdrho));

  tt=supertit(ah,...
	      sprintf('%s ; %s ; %s ; %s\n %s ; %s ; %s',answ1,answ2,answ3,answ4,answ5,answ6,answ7));
  movev(tt,-6)

  % More cosmetics
  set(ah(2),'yaxisl','r')
  
  % Make the plot
  figna=figdisp([],sprintf('%s_%s',Hx,datum),[],0);
  system(sprintf('degs %s.eps',figna));
  system(sprintf('epstopdf %s.eps',figna));
elseif strcmp(Hx,'demo4')
  defval('Gx',[]);
  datum=Gx;
  defval('datum',date)
  
  % The number of parameters to solve for
  np=5;
  % The number of unique entries in an np*np symmetric matrix
  npp=np*(np+1)/2;

  % You need to be in the proper directory to do this, e.g.
  % /u/fjsimons/MyPapers/OLHEDE/Simulations
  f1=sprintf('thzro_%s',datum);
  f3=sprintf('thhat_%s',datum);
  
  % Then load and display what we have
  fid=fopen(f1,'r');
  th0=fscanf(fid,'%f',np)';
  params=fscanf(fid,'%f',8)';
  sclth0=fscanf(fid,'%f',np)';
  likelyjunk=fscanf(fid,'%f',npp)';
  Fisher=fscanf(fid,'%f',npp)';
  truecov=fscanf(fid,'%f',npp)';
  fclose(fid);
  
  thhats=load(f3);
  
  % Now reassemble the covariance
  fullcov=zeros(np,np);
  fullcov(nonzeros(triu(reshape(1:np^2,np,np)')'))=truecov;
  truecov=[tril(fullcov)'+tril(fullcov)-diag(diag(fullcov))];
  % Rescale the covariance matrix
  sclcov=truecov./[sclth0(:)*sclth0(:)'];
  % Actually, best to scale so the diagonal has unit variance
  sclcov=truecov./[diag(sqrt(truecov))*diag(sqrt(truecov))']

  % How about we plot the observed covariance matrix instead
  obscov=cov(thhats);
  obscov=obscov./[diag(sqrt(obscov))*diag(sqrt(obscov))']

  defval('oneortwo',2)

  switch oneortwo
    case 1
     % Plot
     clf
     fig2print(gcf,'portrait')  
     ah=gca;
     imagesc(sclcov)
     kelicol
     caxis([-1 1])
     hold on
     p(1)=plot([2.5 2.5],[0.5 5.5],'k');
     p(2)=plot([0.5 5.5],[2.5 2.5],'k');
     set(p,'Linew',2)
     hold off
     
     % Labels and such
     labs={'D','f2','s2','nu','rho'};
     set(ah,'xtick',1:np,'xtickl',labs)
     set(ah,'ytick',1:np,'ytickl',labs)
     longticks(ah)
     movev(ah,.05)

     shrink(ah,1.25,1.25)
     axis image
     cb=colorbar('hor');
     shrink(cb,1/1.225,1.5)
     movev(cb,-.175); longticks(cb,2)
     axes(cb)
     xlabel('normalized variance/covariance matrix')
     set(cb,'ydir','rev')
     
     t=supertit(ah,sprintf(...
	 '%s on a %ix%i km grid',...
	 sprintf('%i MLE simulations',length(thhats(:,1))),...
	 params(7:8).*params(5:6)/1e3));
     movev(t,.2)
   case 2
     clf
     fig2print(gcf,'portrait')  
     ah=krijetem(subnum(2,2));
     delete(ah(3:4))
     ah=ah(1:2);

     axes(ah(1))
     imagesc(sclcov)
     kelicol
     caxis([-1 1])
     hold on
     p(1)=plot([2.5 2.5],[0.5 5.5],'k');
     p(2)=plot([0.5 5.5],[2.5 2.5],'k');
     set(p,'Linew',2)
     hold off
     axis image
     t(1)=title('predicted');

     axes(ah(2))
     imagesc(obscov)
     kelicol
     caxis([-1 1])
     hold on
     p(3)=plot([2.5 2.5],[0.5 5.5],'k');
     p(4)=plot([0.5 5.5],[2.5 2.5],'k');
     set(p,'Linew',2)
     hold off
     axis image
     t(2)=title('observed');

     set(ah(2),'yaxisl','r')
     movev(ah,-.25)
     serre(ah,1,'across')
     
     % Labels and such
     labs={'D','f2','s2','nu','rho'};
     set(ah,'xtick',1:np,'xtickl',labs)
     set(ah,'ytick',1:np,'ytickl',labs)
     longticks(ah)
     movev(ah,.05)
     
     defval('cblox','ver')
     switch cblox
       case 'hor'
	cb=colorbarf('hor',get(gca,'FontSize'),get(gca,'FontName'),...
		     [0.3 0.3 0.3 0.02]);
	axes(cb)
	xlabel('normalized variance/covariance matrix')
	moveh(cb,0.125)
	movev(cb,-0.02)
      case 'ver'
       cb=colorbarf('ver',get(gca,'FontSize'),get(gca,'FontName'),...
		      [0.175 0.385 0.02 getpos(ah(1),4)]);
       axes(cb)
       ylabel('normalized covariance matrix')
       set(cb,'yaxisloc','l')
     end

     longticks(cb,2)
     
     t=supertit(ah,sprintf(...
	 '%s on a %ix%i km grid',...
	 sprintf('%i MLE simulations',length(thhats(:,1))),...
	 params(7:8).*params(5:6)/1e3));
     movev(t,.5)
  end

  % Make the plot
  figna=figdisp([],sprintf('%s_%s',Hx,datum),[],0);
  system(sprintf('epstopdf %s.eps',figna)); 
end