clear;

%Open binary file for reading ('r') with the precision (i.e. little-endian) of the current system ('n')
fid = fopen('dbz_88d_5km-nmiss.dat','r','n');
testr = fread(fid,'float');
new = reshape(testr,150,150,18);
missing = new == -9999;         %if inequality is true, missing is 1; false, missing is 0.
%new(missing) = NaN;             %assign NaN to places in original array where missing is 1, otherwise use original data
new(missing) = 0.0;             %assign NaN to places in original array where missing is 1, otherwise use original data
dbz = reshape(new,150*150,18);	%[space,time] array 

scale = [-74.5:1.0:74.5];
%contourf(new(:,:,1)');
%figure(2);contourf(new(:,:,2)');

%dimensions for "data" array are [space,time] or [rows,cols]
nspace=size(dbz,1);
ntime=size(dbz,2);

%mininum of 7 points in time domain needed to apply Hilbert (for sin wave problem)
x=[1:nspace];
nx=x;
%t=[0,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25];
t=[28,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50];
nt = [31];

%Plot linear interpolation
data = dbz';
[x,t]=meshgrid(x,t);
[nx,nt]=meshgrid(nx,nt);
inter = interp2(x,t,data,nx,nt,'linear');
linear=reshape(inter,150,150);
%figure(3);contourf(reshape(inter,150,150)')
%contourf(scale,scale,reshape(inter,150,150)')
%set(gca,'clim',[0 50]);
%colorbar



tto=[31:31];
ttn=[28,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50];
%ttn=[0,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25];

sp=nspace;

%Above removes spatial mean, better to do it this way...
[m,n] = size(dbz);
A = mean(dbz,2);	%temporal mean
xc = dbz - repmat(A,1,n); %temporal anomalies for each row
%this method didn't seem to work...
% Apply the Hilbert transform.

test=1;  % test ==1: apply hilber transform on each time series,
         % otherwise, the whole dataset.

if test==1
  data1=[];
  for ii=1:sp
    data1=[data1;hilbert(xc(ii,:))];
  end

else
%this hilbert did work...
  data1=hilbert(xc);
end

%the variance-covariance matrix...where nspace is number of spatial points
covm=data1'*data1/(nspace-1);

[tmp Val Vec]=svd(covm);

Vec2=Vec'*data1';

% Interpolate each time series and reconstruct the data;
%  aa is amplitude ag is phase
recon1=0;

% choose number of modes
nmodes=ntime;
%sum(per(1:nmodes))

for ii=1:nmodes
  tmpt=Vec(:,ii);
  ag=unwrap(angle(tmpt));
  aa=abs(tmpt);
  ag2=interp1(ttn',ag,tto');
  aa2=interp1(ttn',aa,tto');
  tmp=aa2.*exp(i*ag2);
  recon1=recon1+tmp*Vec2(ii,:);
end

recon1=real(recon1);
recon1=recon1';

xc4 = recon1 + A; %temporal anomalies for each row

ceof = reshape(xc4,150,150);

%xc4=recon1;

%contourf(scale,scale,ceof')
%set(gca,'clim',[0 50]);
%colorbar


fid = fopen('goaldbz-step31.dat','r','n');
testrr = fread(fid,'float');
newr = reshape(testrr,150,150);
missing = newr == -9999;         %if inequality is true, missing is 1; false, missing is 0.
%new(missing) = NaN;             %assign NaN to places in original array where missing is 1, otherwise use original data
newr(missing) = 0.0;             %assign NaN to places in original array where missing is 1, otherwise use original data
goal=reshape(newr,150,150);


rmse1 = sqrt(mean(mean((goal-linear).^2)))
rmse2 = sqrt(mean(mean((goal-ceof).^2)))
bias1 = mean(mean(goal-linear))
bias2 = mean(mean(goal-ceof))
