N=100;       % number of time to run the model
K=0.04;      % diffusion coefficient
I=30; J=30;  % domain size

% Velocity field is given below (DO NOT CHANGE)
u=zeros(I,J,1);
v=zeros(I,J,1);
[y_coord,x_coord]=meshgrid(1:I,1:J);
u_adv=1.6;
u=u_adv*cos(3*y_coord/J).*x_coord/(I/2)/1.5;
v=-1.2*u_adv*cos(2*14*x_coord/J);

% how to run the model
% assume you have two sources a t=0
To=zeros(I,J,1);
To(21:22,11:12,1)=5;
To(11:12,15:16,1)=5;
T = adv_diff_model(To,K,N,u,v);

% how to convert between gridded and array
T100 = T(:,:,100);
T_array = T100(:);  % how to convert from gridded to array
T100 = reshape(T_array, [30 30]); % how to convert from array to gridded

% how to make a plot of the model output
subplot(1,2,1); pcolor(x_coord, y_coord, To); shading flat; colorbar;
title 'DISPERSION PATTERN at t=0'
subplot(1,2,2); pcolor(x_coord, y_coord, T100); shading flat; colorbar;
title 'DISPERSION PATTERN at t=100'


%==========================================================
%	how to build E,  for the model y = E*x + n
% NOTE: becuase you know that sources occupy a square of 
% 4 gridpoints, we will build E by perturbing the model 
% initial conditions with impulses over 4 grid points
% square only.
%==========================================================
[I,J]=size(x_coord);
%load E
if ~exist('E')
   n=0;
   E = zeros(I*J, I*J/4);
   for i=1:2:I-1
      for j=1:2:J-1
         n=n+1;
         To=zeros(I,J,1);
         % perturb an inital patch of 4 gridpoints
         To(i:i+1,j:j+1)=1; 
         % run model
         T=adv_diff_model(To,K,N,u,v);
         % get time = 100
         T100=T(:,:,100);
         % save the column of E associated with the patch
         E(:,n)=T100(:);
      end
   end
end

% once you find your estimates of the model parameters x
% you can remap these on the model intial conditions with
% the reverse loop
n=0; To_hat=zeros(I,J,1);
for i=1:2:I-1
   for j=1:2:J-1
      n=n+1;
      To_hat(i:i+1,j:j+1)=xhat(n);
   end
end

% The model
% y=E*yo

% Define a generic cost function 
% J = G(y)

% Infinitesiaml changes in J are given by:
% dJ = dy'*(dG/dy) = (E*dyo)'*(dG/dy) = dyo'*E'*(dG/dy)

% The sensitivity of the cost function to changes 
% in initial conditions dyo becomes:
% dJ/dyo = E'*(dG/dy)
% 
% Going bacl to the class examples we defined 
% J = y'*W*y = G(y)

% where W, the weigths, are choosen so that I am interested in studying 
% the variance only at location 360
[r,c]=size(E);
w=zeros(r,1); w(110)=1; W=w*w';

% Therefore 
% dJ/dyo = E'*(dG/dy) = E'*(2*y'*W)' = 2*E'*W*y

y=ones(r,1);
Jcost=y'*W*y
dJ_dyo=2*E'*W*y;
dJ_dx=2*E'*W*y;


n=0; dJ_dyo_gr=zeros(I,J,1);
for i=1:2:I-1
   for j=1:2:J-1
      n=n+1;
      dJ_dyo_gr(i:i+1,j:j+1)=dJ_dyo(n);
   end
end

% how to make a plot of the model output
subplot(1,2,1); pcolor(x_coord, y_coord, dJ_dyo_gr); shading flat; colorbar;
title 'dJ/dyo'
subplot(1,2,2); pcolor(x_coord, y_coord, reshape(w, [30 30]) ); shading flat; colorbar;
title 'Location where we want to study changes in variance'

% Now assume you have 1 ppm of passive tracer standard deviation at your initial 
% condition and you want to maximize the change at location 110 in the
% future
% so you want std(dJ_dyo*dyo)= 1 =std(dJ_dyo)*dyo
dyo=1/std(dJ_dyo);
yo=dJ_dyo*dyo;

y=E*yo;
Jcost=y'*W*y

n=0; yo_gr=zeros(I,J,1);
for i=1:2:I-1
   for j=1:2:J-1
      n=n+1;
      yo_gr(i:i+1,j:j+1)=yo(n);
   end
end

subplot(1,2,1); pcolor(x_coord, y_coord, yo_gr); shading flat; colorbar;
title 'xhat'
subplot(1,2,2); pcolor(x_coord, y_coord, reshape(y.^2 + 0*w*Jcost, [30 30]) ); shading flat; colorbar;
title 'Value of J at the location'

 

% Now let us verify that this is the best pattern to maximize a change in dJcost by producing
% 1000 random perturbation of the initial condition with a total ammount of passive tracer of 5

for i=1:1000
patt=(randn(c,1)); yo=1/std(patt)*patt;
y=E*yo;
Jcost_ens(i)=y'*W*y;
end
plot(Jcost_ens);
hold on
plot([1 1000], [Jcost Jcost],'k');
legend 'random attemprt to maximize J' 'J maximized using Adjoint'


% Now let us try to just have a point source
for i=1:c
patt=zeros(c,1); patt(i)=1; yo=1/std(patt)*patt;
y=E*yo;
Jcost_single(i)=y'*W*y;
end

clf
plot(Jcost_ens);
hold on
plot(Jcost_single, 'r');
plot([1 1000], [Jcost Jcost],'k');
legend 'random attempt to maximize J' 'Single impulses' 'J maximized using Adjoint'