function Thursday_project;


%
% Decision is to flood a reservoir domain from West or from South. 
% Value is connected to flow volumes. Here we use
% time of flight for value.
%
% Data is processing data to get the variability along north-south
% relative to the variability along east-west. 
%
% Model for permeability field is built from a random parameter beta
% that determines the trend of having a channel zone in the center N-S
% domains. Sometimes this is large, and western flooding is better
% (faster). In other realizations the parameter is small, and southern
% flooding is better (faster).
%
%
% Simulation regression is used to get VOI.
% Realizations of 'permeability gives time of flight, and also the data.
% Simple linear regression is used for E(v|y).
%

n1=25;
n2=25;
n=n1*n2;
% Locations on a grid
sites1=reshape([1:n1]'*ones(1,n2),n,1);
sites2=reshape(ones(n1,1)*[1:n2],n,1);
sites=[sites1 sites2];

figure(1);
clf;
plot(sites1,sites2,'.');

% covariance
a=sites';
c1 = bsxfun(@minus, a(1,:), a(1,:)');
c2 = bsxfun(@minus, a(2,:), a(2,:)');
H=sqrt(c1.^2 + c2.^2);
rr=20; % range
Sigma= 0.3^2*exp(-(3/rr)*H); % exponential

L=chol(Sigma)';

% Monte Carlo sampling

B=1000;
for b=1:B,
    %Prior for channel in middle!

    % Prior in mean
    beta=-gamrnd(1,1,1,1);
    HH=(sites1-n1/2).^2;
    HH=HH/max(max(HH));
    x=HH*beta+L*randn(n,1);
    xm=reshape(x,25,25);
%     figure(10)
%     clf;
%     imagesc(xm);
%     colorbar;
    perm=exp(x);
    perm_m=reshape(perm,n1,n2);
    %figure(2);
    %clf;
    %imagesc(perm_m);
    %colorbar;

    beta
    % Value : Calculate path for time of flight
    
    % Alt1 _ Flooding from the west
    ind1=[13:25:625];
    %perml1=perm(ind1,1); figure(3); plot(perml1)
    tof1=sum(1./perm(ind1));
    val1(b,1)=tof1;

    % Alt2 _ Flooding from the south
    ind2=[12*25+1:26:625];
    %perml2=perm(ind2,1); figure(4); plot(perml2)
    tof2=sum(1.5./perm(ind2));
    val2(b,1)=tof2;
    
    % Data : Calculate ratio of AR(1) correlatioin coefficient along center
    % line N-S and E-W
    vv1=var(perm_m(13,1:25));
    vv2=var(perm_m(1:15,13));
    y(b,1)=log(vv2/vv1);
    disp(sprintf('regr par %0.3g, tof a1 %0.3g, tof a2 %0.3g, var ratio data %0.3g',[beta tof1 tof2 log(vv2/vv1)]));

    hh=1;
end


figure(21);
clf;
plot(y,val1,'kx');
hold;
plot(y,val2,'bo');
hvec=[ones(B,1) y];
[b1]=regress(val1,hvec);
Efit1=hvec*b1;
plot(y,Efit1,'k')
[b2]=regress(val2,hvec);
Efit2=hvec*b2;
plot(y,Efit2,'b')
xlabel('Log ratio of NS , EW variances');
ylabel('Time-of-flight values');


% VOI is PoV- PV
Ev1=mean(val1);
Ev2=mean(val2);
PV=max(-Ev1,-Ev2)

PoV=mean(max([-Efit1 -Efit2]'))

VOI=PoV-PV
