% eda07_04
% example of estimating a filter r that predicts river
% discharge d from precipitation data g; that is, g*r=d
clear all;
global g H;

% in this example, the filter length is chosen to be
% M<<N, to embody the idea that the response of the river
% (meaning the length M of the filter r) is short

M=30; % length of filter r

% load merged discharge and precipitation data
% column 1: time in inter count of days starting 1/1/2002
% column 2: Hudson River discharge in cfs
% column 3: Albany NY precipitation in inches
D = load('precip_discharge_data.txt');

% break out into separate variables, metric conversion
% time in days
t = D(:,1);
% discharge in cubic meters per second
d = D(:,2)/35.3146;
% precipitation in millimeters per day
g = D(:,3)*25.4;

% in this example, we use only a 101 day portion of the
% dataset, mainly so that when we plot it, storm events
% are clearly visible

N=101;
istart = 900;
t = t(istart:istart+N-1);
d = d(istart:istart+N-1);
g = g(istart:istart+N-1);
tmin = min(t);
tmax = max(t);

% the discharge and precipitation data
figure(1);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',3);
axis( [tmin, tmax, -20, 60 ] );
plot( t, g, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('precipitation (mm/day)');
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1000, 3000 ] );
plot( t, d, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('discharge (m3/s)');

% Given the filter equation g*r=d, solve for r given d and g
% using genralized least squares

% Note that the convolution relationship g*r=d is equivalent
% to the matrix equation G r = d where G is the (g*) operator,
% Since we are using filterfun(), we do not need to construct
% G.  We do, however, need to implement the prior information
% equation Hr=h.

% This code can toggle between the prior information of r
% being small, its first derivative being small and its second
% derivative being small, depending on whether DERIVATIVE is 0,
% 1 or 2. (However, the results turn out not to be especially
% sensitive to the choice).
DERIVATIVE=1;

if( DERIVATIVE==2 )
% The prior information equation H*r = h is implemented
% for the smoothness condition of the second derivative
% being close to zero
e=2*max(abs(g)); % the coefficient was manually tuned
                   % to produce a reasonably smooth filter
K=M;
L=N+K;
H=spalloc(K,M,3*K);
for j = [2:K]
    H(j,j-1)=e;
end
for j = [1:K]
    H(j,j)=-2*e;
end
for j = [1:K-1]
    H(j,j+1)=e;
end
h=zeros(K,1);

elseif (DERIVATIVE==1)

% The prior information equation H*r = h is implemented
% for the smoothness condition of the first derivative
% being close to zero
e=4.0*max(abs(g)); % the coefficient was manually tuned
                   % to produce a reasonably smooth filter
K=M;
L=N+K;
H=spalloc(K,M,3*K);
for j = [1:K]
    H(j,j)=e;
end
for j = [1:K-1]
    H(j,j+1)=-e;
end
h=zeros(K,1);

else % DERRIVATIVE == 0
    
% The prior information equation H*r = h is implemented
% for the smallness condition of the filter is close to zero
e=0.1*max(abs(g)); % the coefficient was manually tuned
                   % to produce a reasonably small filter
K=M;
L=N+K;
H=spalloc(K,M,3*K);
for j = [1:K]
    H(j,j)=e;
end
h=zeros(K,1);
end

% The generalized least squares equation is Fr=f
% with F=[G;H] and f=[d;h]
f(1:N)=d;
f(N+1:N+K)=h;

% set up the right hand side F'f = GT d + HT h
% note that GT v is similar to the convolution G v = g*v
% except that p is time-reversed and position of
% results in convolution is shifted to the bottom
temp=conv(flipud(g),d);
FTfa = temp(N:N+M-1);
FTfb = H'*h;
FTf=FTfa+FTfb;

% solve
r=bicg(@filterfun,FTf,1e-10,3*L);

% plot the filter r
figure(2);
clf;
set(gca,'LineWidth',3);
hold on;
axis( [0, M, -5, 10 ] );
plot( [1:M]', r, 'k-', 'LineWidth', 3 );
plot( [1,M]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('filter r');

% plot the predicted discharge
dpre = conv(g,r);
dpre = dpre(1:N);
figure(1);
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1000, 3000 ] );
plot( t, dpre, 'k-', 'LineWidth', 2 );





