spatialPatternAR1(DIM,PHI,SIGMA,PERIODIC,GENERATOR),

function x = spatialPatternAR1(DIM,PHI,SIGMA,PERIODIC,GENERATOR),
% function x = spatialPatternAR1(DIM,PHI,SIGMA,PERIODIC,GENERATOR),
%
% This function generates a first-order autoregressive spatial process with
% a normal error distribution
%
%     DIM is a two component vector that sets the size of the spatial pattern
%            The routine only works for spatial data, so all componentes of
%            DIM must be gretaer than 1.
%           (DIM=[10,5] is a 10x5 spatial grid)
%     PHI is the variable governing the degree of autocorrelation 
%           (PHI must be less than 0.25)
%            PHI = 0    No autocorrelation
%            PHI = 0.25 A uniform plane with no error possible.
%     SIGMA is the vaurance of the normally distributed error
%           Default is a variance of 1
%     PERIODIC = 1 then the spatial pattern is periodic
%              = 0 then the spatial pattern has fixed boundaries
%            Default is PERIODIC = 1
%     GENERATOR = 0 the return a realisation of the spatial process
%               = 1 then return a matrix which can be used to quickly
%                   generate an AR1 spatial process
%            Default is GENERATOR = 0
%            If many spatial patterns are being generated, it is much
%            faster to calculate the generator matrix once only.

% The autoregressive process is  
%      X_i = PHI * (X(i-1,j) + X(i+1,j) +X(i,j-1) +X(i,j+1)) + error
%
% the matrix M is inv(I-PHI*W) where W is the first order
% neighbourhood matrices
%                            X(i-1,j)
%                                |
%                  X(i,j-1) -  X(i,j) - X(i,j+1)
%                                |
%                             X(i+1,j)

% To generate autocorrelated data,
%  1/ Start with a independent random variable for each grid location with a
%       given error distribution, e.
%  2/ Rearrange this into a column vector, e_vec=reshape(e,prod(DIM),1)
%  3/ If x_vec is a random variable with the required autoregressive property.
%     Then we must solve the equation
%                M * x_vec = e_vec
%      by using Gaussian quadrature (this is computationally efficient)
%      so that
%             x_vec = M \ e_vec
%     The inverse of matrix M, is the matrix output by the function when
%     GENERATOR = 1. In this case the autoregressive data is generated by 
%     x_vec = inv(M) * e_vec
%  4/ To regain the data in grid form, x= reshape(x_vec,dim)
% Note that the error distribution is not restricted to be normal. e_vec
% can have any desired error distribution.

% Everything is calculated using sparse matrices and then converted to full
% matrices at the end

% Written by Jon Yearsley 1 may 2004
%            [e-mail]

if nargin<2,
    error('Function requires at least 2 input arguments')
elseif nargin <3,
    SIGMA = 1;
    PERIODIC = 1;
    GENERATOR = 0;
elseif nargin<3,
    % Assume cyclical boundary conditions
    PERIODIC = 1;
    GENERATOR = 0;
elseif nargin<4,
    GENERATOR = 0;
end

if PHI>=0.25,
    error('PHI must be less than 0.25')
    return;
end

if nnz(DIM<=1) | length(DIM)~=2
    error('DIM must be a two component vector, and all components must the greater than 1.')
    return;
end

if SIGMA<0,
    error('SIGMA must not be negative')
    return;
end

N = prod(DIM);
W = sparse(N,N);

if PERIODIC,
    template = sparse(DIM(1),DIM(2));
    template(1,2)=1;
    template(2,1)=1;
    template(1,DIM(2))=1;
    template(DIM(1),1)=1;
end

% Create the first-order neighbourhood matrix.
% note that any order of autoregressive process can be simulated using this
% procedure if the correct neighbourhood matrices are created
for i=1:N,
    if PERIODIC,
        shift = [mod(i-1,DIM(1)) , floor((i-1)/DIM(1))];
        W(i,:) = reshape(circshift(template,shift),1,N);
    else
        % Check we're not at the start of a row
        if mod(i,DIM(1))~=1, W(i,i-1) = 1; end;
        % Check we're not at the end of a row
        if mod(i,DIM(1))~=0, W(i,i+1) = 1; end;
        % Check we're not at the start of a column
        if floor((i-1)/DIM(1))>0, W(i,i-DIM(1)) = 1; end;
        % Check we're not at the end of a column
        if floor((i-1)/DIM(1))<DIM(1)-1, W(i,i+DIM(1)) = 1; end
    end
end

% Create the generator matrix
M = speye(N) - PHI*W;

% Now generate the spatial pattern with normal error structure
if GENERATOR,
    x = inv(M);
else
    x = reshape(M\randn(N,1)*sqrt(SIGMA),DIM);
end