Cokriging Weights

1. template <
             class MatrixLibrary,
             class Location,
             class InputIterator,
             class CovarianceSet,
             class KrigingConstraints,
             class Vector
            >
   int
   cokriging_weights(
                     Vector& weights, double kriging_variance, 
                     const Location& center,
                     InputIterator first_neigh, InputIterator last_neigh,
                     CovarianceSet& covar,
                     KrigingConstraints& Kconstraints
                     ) 
   
   
   
2. template <
             class MatrixLibrary,
             class Location,
             class InputIterator,
             class CovarianceSet,
             class KrigingConstraints,
             class Vector
            >
   int
   cokriging_weights(
                     Vector& weights, 
                     const Location& center,
                     InputIterator first_neigh, InputIterator last_neigh,
                     CovarianceSet& covar,
                     KrigingConstraints& Kconstraints
                     ) 
   
   
   

The cokriging_weights algorithm solves a cokriging system :

$$\left\{\vphantom{ \begin{array}{l} Var \Big( \sum_{\alpha=1}^{n}... ..._{\alpha_1} \},\ldots,\{\lambda_{\alpha_{N_v}} \} \Big)=0 \end{array} }\right.$$$$\begin{array}{l} Var \Big( \sum_{\alpha=1}^{n} \lambda_{\alpha} ... ...\{\lambda_{\alpha_1} \},\ldots,\{\lambda_{\alpha_{N_v}} \} \Big)=0 \end{array}$$

where f_i$\Big(${$\lambda_{\alpha}^{}$},{$\lambda_{\alpha_1}^{}$},...,{$\lambda_{\alpha_{N_v}}^{}$}$\Big)$, i = 1,..., p are p constraints expressed by the Kconstraint. The solution to this system is a set of weights:

($\lambda_{1}^{}$,...,$\lambda_{n}^{}$,$\lambda_{1}^{1}$,...,$\lambda^{N_v}_{n_{\tiny N_v}}$,$\mu_{1}^{}$,...,$\mu_{p}^{}$), where the weights $\mu_{j}^{}$ are the Lagrange weights used to account for the constraints f₁,..., f_p.

The cokriging_weights function stores the kriging weights into Vector weights in the following order:

($\lambda_{1}^{}$,...,$\lambda_{n}^{}$,$\lambda_{1}^{1}$,...,$\lambda^{N_v}_{n_{\tiny N_v}}$,$\mu_{1}^{}$,...,$\mu_{p}^{}$).

Only the weights $\lambda_{1}^{}$,...,$\lambda^{N_v}_{n_{\tiny N_v}}$ are useful to compute the kriging estimate. The weights $\mu_{1}^{}$,...,$\mu_{p}^{}$ are used to compute the kriging variance. Vector weights does not need to be of the correct size n + p when passed to the function.

Two versions of the algorithm exist. Version 1 computes the kriging variance and stores it into kriging_variance, while Version 2 does not compute the kriging variance (there are cases where the kriging variance is useless, e.g. in indicator kriging). Both version of the algorithm allow to change the linear algebra library (which defines matrices, matrix inversion routines, ...) used by cokriging_weights (see 2.4). The default linear algebra library is the TNTlibrary (slightly modified), a public domain library by Roldan Pozo, Mathematical and Computational Sciences Division,National Institute of Standards and Technology.

The value returned by kriging_weights indicates if any problem were encountered:

• returns 0 if the function executed successfully
• returns 1 if the kriging system could not be solved
• returns 2 if the neighborhood of conditioning data was empty

Where defined

In header file <kriging.h>

Requirements on types

• Vector is a container on which an iterator is defined. It must have three member functions whose prototypes are:
  1. iterator Vector::begin() which returns an iterator to the first element of the Vector
  2. size_type Vector::size() which returns the size of the Vector.
  3. void Vector::resize(size_type n) which changes the size of the Vector to n
• Location is a model of Location.
• InputIterator is a model of Input Iterator (see STL requirements). It iterates on a range of Neighborhood.
• CovarianceSet is a model of Covariance Set. In expression covar(i,j,u1,u2), required by Covariance Set, u1 and u2 must be of type Location.
• KrigingConstraints is a model of Kriging Constraint.
• MatrixLibrary specifies the linear algebra library to use. The requirements on matrix_lib are fully defined in 2.4. By default, it is the TNTlibrary (slightly modified), a public domain library by Roldan Pozo, Mathematical and Computational Sciences Division,National Institute of Standards and Technology.

Remarks

Vector weights is resized after it is passed to cokriging_weights, unless it is of the correct size. The cost of this resize can be decreased by smartly managing the Vector's memory (in the style of an STL vector). However, each call to cokriging_weights implies solving a linear system of equations, hence the cost of the resize would usually be negligible.

Example

Estimate z(0, 0) from z(2, 3) = 0.21, z(4, - 7) = 0.09 and one secondary data s(1, 3) = 42.1 by ordinary (full-)cokriging. After the call to kriging_weight, vector weights has size 5; its three first elements are the weights corresponding to z(2, 3), z(4, - 7) and s(1, 3) respectively, and its two last elements are the Lagrange parameters.

typedef std::vector<geo_value2d> neighborhood;

int main()
{  
  location2d u1(2,3);
  location2d u2(4,-7);

  geo_value2d Z1(u1,0.21);
  geo_value2d Z2(u2,0.09); 
  
  Neighborhood neighbors1, neighbors2;
  neighbors1.push_back(Z1);
  neighbors1.push_back(Z2);

  location2d u3(1,2);
  geo_value2d S1(u3,42.1);
  neighbors2.push_back(S1);

  neighborhood* neigh_array[2]={&neighbors1, &neighbors2};

  location2d u(0,0);

  std::vector<double> weights;

  double kvariance;
  kriging_weights(weights, kvariance,
                  u, neigh_array, 2,
                  LMC_covariance, OK_constraint);
}