function Friday_project;


% Forestry decision, harvest all cells ( 25 x 25 on a grid) or none. 
% Data along survey lines help make decision, but data has a price. 
% Sequential information gathering at N-S lines (length 25).
% Myopic search of N-S lines in the 25 x 25 grid.
%
% Prior is a Gaussian process.
% Data is unbiased with Gaussian noise.
%
% File shows the results of the myopic strategy, for one realization
% of data. The results will vary from realization to realization. 
% Sometimes many lines are surveyed, sometimes only a few. 
% If posterior value is desired this is obtained by averaging over several such runs.
%

% 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;
mx=m*ones(n,1);

% 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;
Sigma= r^2*exp(-(3/range)*H); % exponential
% Cholesky factor
L=chol(Sigma)';
figure(2)
clf;
imagesc(Sigma);

% Likelihood noise
tau=10;
% dim of measurement
M=25;

% Price of experiment
Price=0.1;

strule=0;
cnt=1;
while (strule==0)
    
    if (cnt==1)
        Stop=0;
    else
       Stop=max(sum(mx),0); 
    end
    disp(sprintf('Stop value is %0.3g.',Stop));

    Cont=zeros(n1,1);
    maxCont=-10000;
    maxind=0;
    for i=1:n1,
        % Set design points
        F=zeros(M,n);
        F(1:M,(i-1)*n1+1:i*n1)=eye(M);
 
        R=Sigma*F'*((F*Sigma*F'+tau^2*eye(M))\F*Sigma);

        mw=sum(mx);
        rw=sqrt(sum(sum(R)));
        aw=mw/rw;
        if (cnt==1)
            % VOI, no comp with price. 
            Cont(i)=mw*normcdf(aw)+rw*normpdf(aw)-max(0,mw);
        else
            % Additional PostValue, at the price of more data
            Cont(i)=mw*normcdf(aw)+rw*normpdf(aw)-Price;
     
        end
        if (Cont(i)>maxCont)
            maxCont=Cont(i);
            maxind=i;
        end
        %disp(sprintf('Cont value line %d is %0.3g.',i,Cont(i)));

    end

    if (maxCont<Stop)
        strule=1;
        disp('Stop');
    else
        disp(sprintf('Continue to step %d. Sampled line %d',cnt+1,maxind));
        F=zeros(M,n);
        F(1:25,(maxind-1)*n1+1:maxind*n1)=eye(25);
 
        C=F*Sigma*F'+tau^2*eye(M);

        % sample data
        y=F*mx+chol(C)'*randn(M,1);
        
        % Update model
        xhat=mx+Sigma*F'*(C\(y-F*mx));
        mx=xhat;
        % Kriging variance
        vhat=Sigma-Sigma*F'*(C\(F*Sigma));
        Sigma=vhat;
        figure(1);
        clf;
        imagesc(reshape(diag(vhat),n1,n1)); colorbar;
        pause(1);
    end
    
    cnt=cnt+1;
end

