function [x0s,s,sv,id,l]=cokri(x,x0,model,c,itype,avg,block,nd,ival,...
                               nk,rad,ntok,b)
%
% COKRI performs point or block cokriging in D dimensions (any integer)
%       of P variables (any integer) with a combination of K basic models
%       (any integer).
%
% Syntax:
% [x0s,s,sv,id,l]=cokri(x,x0,model,c,itype,avg,block,nd,ival,nk,rad,ntok)
%
% Input description:
% ------------------
%   x: The n x (p+d) data matrix. This data matrix can be imported from an
%      existing ascii file. Missing values are coded 'nan' (not-a-number).
%   x0: The m x d matrix of coordinates of points to estimate.
%   model: Each row of this matrix describes a different elementary structure.
%          The first column is a code for the model type, the d following
%          columns give the ranges along the different coordinates and the
%          subsequent columns give rotation angles (a maximum of three).
%          The codes for the current models are:
%             1: nugget effect
%             2: exponential model
%             3: gaussian model
%             4: spherical model
%             5: linear model
%             6: quadratic
%             7: power (h^b)
%             8: logarithmic
%             9: sinc; i.e. sin(h) / h
%            10: Bessel [ Jo (h) ]
%            11: exp(-h)   * cos(bh)
%            12: exp(-h)   *  Jo(bh)
%            13: exp(-h^2) * cos(bh)
%            14: exp(-h^2) *  Jo(bh)
%	     15: exp(-h^2) * (1 - br^2)
%
%   *** a new parameter, b, has been added to cokri and cokri2 ***
%
%          Note: a linear model is specified by arbitrary ranges and a sill
%                such that sill/range gives the desired slope in the direction
%                condidered.
%   c: The (rp x p) coefficient matrix of the coregionalization model.
%      Position (i,j) in each submatrix of size p x p give the sill of the
%      elementary component for each cross-variogram (variogram) between
%      variable i and variable j.
%   itype: Code to indicate which type of cokriging is to be performed:
%             1:  simple cokriging
%             2:  ordinary cokriging with one nonbias condition
%                 (Isaaks and Srivastava).
%             3:  ordinary cokriging with p nonbias condition.
%             4:  universal cokriging with drift of order 1.
%             5:  universal cokriging with drift of order 2.
%             99: cokriging is not performed, only sv is computed.
%   block: Vector (1 x d), giving the size of the block to estimate;
%          any values when point cokriging is required.
%   nd: Vector (1 x d), giving the discretization grid for block cokriging;
%       put every element equal to 1 for point cokriging.
%   ival: Code for cross-validation.
%             0:  no cross-validation
%             1:  cross-validation is performed  by removing one variable at a
%                 time at a given location.
%             2:  cross-validation is performed by removing all variables at a
%                 given location.
%   nk: Number of nearest neighbors in x matrix to use in the cokriging
%       (this includes locations with missing values even if all variables
%       are missing).
%   rad: Search radius for neighbors.
%   ntok: Points in x0 will be kriged by groups of ntok grid points.
%         When ntok>1, the search will find the nk nearest samples within
%         distance rad from the current ntok grid points centroid.
%컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴컴
% Output description:
% -------------------
%
%   For the usual application, only x0s and s are required and the other
%   output matrices may be omitted.
%
%   x0s: m x (d+p) matrix of the m points (blocks) to estimate by the
%        d coordinates and p cokriged estimates.
%   s: m x (d+p) matrix of the m points (blocks) to estimate by the
%      d coordinates and the p cokriging variances.
%   sv: 1 x p vector of variances of points (blocks) in the universe.
%   id: (nk x p) x 2 matrix giving the identifiers of the lambda weights for
%       the last cokriging system solved.
%   l: ((nk x p) + nc) x (ntok x p) matrix with lambda weights and
%      Lagrange multipliers of the last cokriging system solved.

%
casesen off;

% definition of some constants.

[m,d]=size(x0);

% check for cross-validation

if ival>=1,
   ntok=1;
   x0=x(:,1:d);
   nd=ones(1,d);
   [m,d]=size(x0);
end
[rp,p]=size(c);
[n,t]=size(x);
nk=min(nk,n);
ntok=min(ntok,m);
idp=[1:p]';
ng=prod(nd);

% echo of data and input options and parameters

%echo_input=['dimension ',num2str(d)]
%echo_input=['number of variables ',num2str(p)]
%'Strike any key when ready',pause
%x
%['Strike any key when ready'],pause
%x0
%['Strike any key when ready'],pause
%model
%c
%['Strike any key when ready'],pause
%itype,ival
%avg
%block,nd
%['Strike any key when ready'],pause
%nk,rad,ntok

% compute point (ng=1) or block (ng>1) variance

for i=1:d,
   nl=prod(nd(1:i-1));
   nr=prod(nd(i+1:d));
   t=[.5*(1/nd(i)-1):1/nd(i):.5*(1-1/nd(i))]';
   t2=[t2,kron(ones(nl,1),kron(t,ones(nr,1)))];
end
grid=t2.*(ones(ng,1)*block);
t=[grid,zeros(ng,p)];

% for block cokriging, a double grid is created by shifting slightly the
% original grid to avoid the zero distance effect (Journel and Huijbregts, p.96).

if ng>1,
   grid=grid+ones(ng,1)*block/(ng*1e6);
end
[x0s,s,id,l,k0]=cokri2(t,grid,[],model,c,sv,99,avg,ng,b);

% sv contain the variance of points or blocks in the universe

for i=1:p,
   sv=[sv,means(means(k0(i:p:ng*p,i:p:ng*p))')];
end

% start cokriging

for i=1:ntok:m,
   nnx=min(m-i+1,ntok);
%   ['kriging points #',num2str(i),' to ', num2str(i+nnx-1)]

   % sort x samples in increasing distance relatively to centroid of 'ntok'
   % points to krige

   centx0=ones(n,1)*means(x0(i:i+nnx-1,:));
   tx=[x(:,1:d)-centx0].*[x(:,1:d)-centx0]*ones(d,1);
   [tx,j]=sort(tx);

   % keep samples inside search radius; create an identifier of each sample
   % and variable (id)

   t=[];
   id=[];
   ii=1;
   tx=[tx;nan];
   while ii<=nk & tx(ii)<rad*rad,
      t=[t;x(j(ii),:)];
      id=[id;[ones(p,1)*j(ii),idp]];
      ii=ii+1;
   end
   t2=x0(i:i+nnx-1,:);

   %  if block cokriging discretize the block

   t2=kron(t2,ones(ng,1))-kron(ones(nnx,1),grid);

   % check for cross-validation

   if ival>=1,
      est=zeros(1,p);
      sest=zeros(1,p);

      % each variable is cokriged in its turn

      if ival==1,
         np=1;
      else
         np=p;
      end
      for ip=1:np:p,

         % because of the sort, the closest sample is the sample to
         % cross-validate and its value is in row 1 of t; a temporary vector
         % keeps the original values before performing cokriging.

         vtemp=t(1,d+ip:d+ip+np-1);
         t(1,d+ip:d+ip+np-1)=ones(1,np)*nan;
         [x0ss,ss]=cokri2(t,t2,id,model,c,sv,itype,avg,ng,b);
         est(ip:ip+np-1)=x0ss(ip:ip+np-1);
         sest(ip:ip+np-1)=ss(ip:ip+np-1);
         t(1,d+ip:d+ip+np-1)=vtemp;
      end
      x0s=[x0s;[t2,est]];
      s=[s;[t2,sest]];
   else
      [x0ss,ss,id,l]=cokri2(t,t2,id,model,c,sv,itype,avg,ng,b);
      x0s=[x0s;[x0(i:i+nnx-1,:),x0ss]];
      s=[s;[x0(i:i+nnx-1,:),ss]];
   end
end

casesen on;
