function NorwegianWood;

%
%  Forestry example, 
% Farmer must decide whether to harvest a forest, or not. 
% Low decision flexibility, all or nothing, no selection of sub-units.
%
% Surveying is done at M random locations.
%
% VOI is analytically available, and we can easily do sensitivity analysis
% of the different model parameters (mean, corr range, etc).
%

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

% Display the grid nodes
figure(1);
clf;
plot(sites1,sites2,'.');

% Prior mean
m=0;

% covariance is computed from distances
a=sites';
c1 = bsxfun(@minus, a(1,:), a(1,:)');
c2 = bsxfun(@minus, a(2,:), a(2,:)');
H=sqrt(c1.^2 + c2.^2);
range=20; % range
r=1;
phi=(3/range);
Sigma= r^2*(1+phi*H).*exp(-phi*H); % exponential
% Cholesky factor
L=chol(Sigma)';
figure(2)
clf;
imagesc(Sigma);

% Sample from the gaussian process
x=L*randn(n,1);
mx=m*ones(n,1);
x=x+mx;
xm=reshape(x,n1,n1);
% Image of sample
figure(3)
clf;
imagesc(xm);
colorbar;
% Set design points
M=10;
des=ceil(n*rand(M,1));
F=zeros(M,n);
for i=1:M,
   F(i,des(i))=1;
end
% Measurement st dev
tau=0.01;
% sample data
y=F*x+tau*randn(M,1);

%figure(4);
%clf;
%plot(y)

% Kriging prediction
xhat=mx+Sigma*F'*((F*Sigma*F'+tau^2*eye(M))\(y-F*mx));
xhatmat=reshape(xhat,n1,n1);
figure(5);
clf;
imagesc(xhatmat); colorbar;

% Kriging variance
vhat=diag(Sigma-Sigma*F'*((F*Sigma*F'+tau^2*eye(M))\(F*Sigma)));
figure(6);
clf;
imagesc(reshape(vhat,n1,n1)); colorbar;

gg=0;

% Reduction in variance is R=Sigma*F'*((F*Sigma*F'+tau^2*eye(M))\F*Sigma);
% The VOI is then using the m/r
% in the cumulative and probability density

R=Sigma*F'*((F*Sigma*F'+tau^2*eye(M))\F*Sigma);

mw=sum(mx);
rw=sqrt(sum(sum(R)));
aw=mw/rw;
VOI=mw*normcdf(aw)+rw*normpdf(aw)-max(0,mw);

% Sensitivity to mu
mvec=[-1:0.1:1];
VOIvec=zeros(1,length(mvec));
for i=1:length(mvec),
   mx=mvec(i)*ones(n,1); 
   mw=sum(mx);
   aw=mw/rw;
    VOIvec(i)=mw*normcdf(aw)+rw*normpdf(aw)-max(0,mw);

end
figure(10);
clf;
plot(mvec,VOIvec)
xlabel('Prior mean level'); ylabel('VOI');

mx=zeros(n,1);

% Sensitivity to corr range
rangevec=[5:5:50];
VOIvec=zeros(1,length(rangevec));
for i=1:length(rangevec),
    % Cov matrix with this range
    range=rangevec(i);
    Sigma= r^2*exp(-(3/range)*H); % exponential
    R=Sigma*F'*((F*Sigma*F'+tau^2*eye(M))\F*Sigma);
    rw=sqrt(sum(sum(R)));
    aw=mw/rw;
    mw=sum(mx);
    aw=mw/rw;
    VOIvec(i)=mw*normcdf(aw)+rw*normpdf(aw)-max(0,mw);

end
figure(11)
clf;
plot(rangevec,VOIvec);
xlabel('Spatial correlation range'); ylabel('VOI');
range=40;
