function [vgrid,vgriderror,vdata,vdataerror] = ...
    multi_oi_err(data,posgrid,posdata,regions,lonscl,latscl,timscl,e)
% Optimal interpolation of 3-d dataset
%
% [vgrid,vgriderror,vdata,vdataerror] = ...
%        multi_oi_err(data,posgrid,posdata,regions,lonscl,latscl,timscl,e)
%
% Inputs:
%
% DATA     m1*1 vector of data
% POSGRID  m2*3 matrix of grid locations [usually latitude,longitude,time]
% POSDATA  m1*3 matrix of data locations [usually latitude,longitude,time]
%
%          (m1 is the number of data points)
%          (m2 is the number of grid points)
%
% REGIONS  n*5 matrix defining mapping regions [lonW,lonE,latS,latN,overlap]
% LONSCL   scalar longitude scale (in degrees)
% LATSCL   scalar latitude scale (in degrees)
% TIMSCL   scalar time scale (in days)
% E        scalar noise to signal variance ratio
% 
% Outputs:
% 
% VGRID       m2*1 vector of mapped fields
% VGRIDERROR  m2*1 estimated error to signal ratio (assumed statistics)
% VDATA       m1*1 vector of mapped fields on original data locations
% VDATAERROR  m1*1 estimated error to signal ratio of the mapped data
%
% If outputs are not requested they are not computed, and this runs a lot
% faster. 
%
%
% Purpose:
%
%    An optimal mapping routine, taking data measured in arbitrary
%    geographic locations and mapping these data onto a more regular
%    grid.  This routine is designed to handle the common problems
%    associated with large data sets, by mapping only in predetermined
%    regions, rather than using all the data globally.  As should 
%    happen in every mapping problem the data are both mapped onto
%    the prescribed grid, and the data locations (so that the mapped
%    field can be checked with the original data to ensure that the 
%    statistics are valid and consistent).
%
% Notes for the apprehensive:
%
%    Regions matrix allows the data set has to be broken into small 
%    regions.  The regions are defined as described in function 
%    GENSUBREGIONS.  The overlap parameter controls the amount of overlap 
%    of each of the data regions.  It will help eliminate edge effects at 
%    the boundaries of the regions provided it is large enough that the 
%    covarianvce function becomes small over the overlap distance.
%    But if the overlap parameter is too large, the number of data will 
%    increase dramatically, slowing the code.  In problems
%    where there are less than 1000 data pts and grid points it is simpler
%    to define just one region that includes all the data.
%
%    The data do not have to be sorted into any order, and the returned data
%    correspond to the ordering specified in the posgrid and posdata 
%    variables.
%
%    Possible improvements would be to allow variable (i) signal, noise,
%    latitude, longitude and overlap scales.
%
% Author: N.L Bindoff
% Date: 24 April 1994
% 
% Mods: A. Walker and J. Wilkin, July 1997

% Set aside the memory 
vdata = zeros([size(posdata,1) 1]);
vdataerror = vdata;
vgrid = zeros([size(posgrid,1) 1]);
vgriderror = vgrid;

% Create the data-data covariance function.  In this case it is the
% same for all data so it is done just once.

totaltime=cputime;

for ii = 1:size(regions,1)
  regiontime = cputime;
 
  gridindex = findinrange(posgrid(:,[2 1]),regions(ii,1:4));
  dataindex = findinrange(posdata(:,[2 1]),...
      regions(ii,1:4) + regions(ii,5)*[-1 1 -1 1]);
  disp(['   subregion ' num2str(ii) ': ' num2str(length(gridindex)) ...
      ' grid points and ' num2str(length(dataindex)) ' data points']);

  if ~isempty(dataindex)

    % find signal standard deviation s from data(dataindex)*data(dataindex)
    s = std(data(dataindex(:)));

    % create data-data covariance matrix
    Cdd = coval_yxt(...
	posdata(dataindex,:),posdata(dataindex,:),lonscl,latscl,timscl);
    if e~=0
      mm = size(dataindex,1);
      lam = e/s;
      Cdd(1:mm+1:mm^2+mm) = diag(Cdd)+lam.*ones([mm 1]);  
    end
    
    % Calculate the objectively mapped fields on the data and grid fields.
    % Note the mean field is removed to ensure that estimate is unbiased
    meanv = mean(data(dataindex));
    wght = Cdd\(data(dataindex)-meanv.*ones([length(dataindex) 1]));
    
    if nargout > 1 % only invert Cdd if error estimates were requested
      Cdd=inv(Cdd);		
    end
    
    % Mapping the data back onto itself
    if nargout > 2 % only map the data if requested 
      for j=1:size(dataindex,1)
	Cmd = coval_yxt(posdata(dataindex(j),:),...
	    posdata(dataindex,:),lonscl,latscl,timscl);
	vdata(dataindex(j)) = Cmd*wght+ meanv;
	if nargout > 3 % only compute the error if requested
	  vdataerror(dataindex(j)) = sqrt(1-Cmd*Cdd*Cmd'); % include a factor
	  % of s inside the sqrt to get a dimensional error estimate
	end
      end
    end
    
    % Mapping onto the model grid
    for j=1:size(gridindex,1)
      [Cmd] = coval_yxt(posgrid(gridindex(j),:),...
	  posdata(dataindex,:),lonscl,latscl,timscl);
      [vgrid(gridindex(j))] = Cmd*wght + meanv;
      if nargout > 1 % only compute the error if requested
	[vgriderror(gridindex(j))] = sqrt(1-Cmd*Cdd*Cmd'); % include a factor
	% of s inside the sqrt to get a dimensional error estimate
      end
    end
  end
end
