function anal = enkf(dat,HHH,vr,fore);
% EKFmem(dat, HHH, vr, fore)
%   A simple Ensemble Kalman Filter
%   Assumes that DAT and VR are the same dimensions
%   as the state vector (build from MEM).
%
% INPUTS
%   DAT   is the data for the update
%   HHH   is the observation matrix (typically 0's and 1's)
%   VR    is a vector of observation variances.
%   FORE  is a matrix of the forecast ensemble, columns are replicates
%
% OUTPUTS
%   Updates the elements in MEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3

Nens = size(fore,2);
Nx   = size(fore,1);
rH   = size(HHH,1);
cH   = size(HHH,2);
dat  = dat(:);
vr   = vr(:);
Ny   = length(dat);
if(rH ~= Ny); error('Dimension of Data is incorrect'); end
if(Ny ~= length(vr)); 
    error('Dimension of Variance vector is incorrect'); 
end


%rrr = perfectnorm(Ny,Nens,1);
rrr = randn(Ny,Nens);
ddd = dat*ones(1,Nens) + (sqrt(vr)*ones(1,Nens)).*rrr;
rrr = ddd - HHH*fore; 

mne = mean(fore')';
eee = fore - mne*ones(1,Nens);
if(Nx < 100)
  CCC = cov(rrr');
  PPP = HHH*CCC*HHH' + diag(vr);
  rrr = CCC*(HHH/PPP)*rrr;
else
  hhe = HHH*eee;
  PPP = hhe*(hhe'./(Nens-1)) + diag(vr);
  rrr = PPP\rrr;
  rrr = hhe'*rrr;
  rrr = (eee./(Nens-1))*rrr;
end
anal = fore + rrr;