% eda07_05
% prediction error filter (pef) for Neuse Rievr Hydrograph
% this code uses generalized least squares to implement
% the prior information that pef(1)=(-1).  Two solutions
% are calculated, one via standard matrix algebra,
% (FTF)\(FTf) and the other via the biconjugate gradient method

clear g H;
global g H;

% read in data
D = load('neuse.txt');

% set up timeseries
t=D(:,1);
d=D(:,2);
N = length(d);
Dt = t(2)-t(1);
tmin = 0.0;
t = tmin + Dt*[0:N-1]';
tmax = t(N);

% sizes of various arrays
M=100; % length of pef is user tunable
K=1;
L=N+K;

% for a prediction error filter, pef, the data is the filter, g
g = d;

% prior information that first element of the pef is (-1)
e=(10^1)*max(abs(g)); % large weight
H=zeros(K,M);
H(1,1)=e;
h=zeros(1,1);
h(1)=(-e);

% G needed only for standard colution
G=zeros(N,M);
G=toeplitz(g, [g(1), zeros(1,M-1)] );

% F needed only for standard solutiom
F=zeros(L,M);
F(1:N,:)=G;
F(N+1:N+K,:)=H;

% f needed only for standard solution
f=zeros(L,1);
f(1:N)=zeros(N,1);
f(N+1:N+K)=h;

% pef by standard matrix algebra
pef1=(F'*F)\(F'*f);

% FT f needed only for biconjugate gradient method
% FTf = GT 0 + HT h = [e 0 0 ... 0]T [-e]
FTf = zeros(M,1);
FTf(1) = -e^2;

% solve for predictione error filter using generalized least squares
% constraint that pef(1)=(-1) implemented with generalized least squares
pef2=bicg(@filterfun,FTf,1e-12,3*L);

% prediction error
perror = conv(pef2,g);

% plot pef
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [ 0, Dt*(M-1), min(pef1), max(pef1)] );
plot( t(1:M), pef1, 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'pef(t) (std)' );
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [ 0, Dt*(M-1), min(pef1), max(pef1)] );
plot( t(1:M), pef2, 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'pef(t)( bicg)' );

% plot hydrograph data and prediction error for first year
figure(2);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
Np = 365;
axis( [tmin, tmin+Dt*Np, -5000, 15000]' );
hold on;
plot( t(1:Np), d(1:Np), 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'd(t)' );
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmin+Dt*Np, -5000, 15000]' );
plot( t(1:Np), perror(1:Np), 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'e(t)' );








