GsTL介绍
1、为什么要有GsTL
简单来讲GsTL（Geostatistal Template Library）就是大名鼎鼎的GSLIB的C++版本，由于GSLIB是fortran编写的，越来越不适应当代编程的需要，因此斯坦福大学的SCRF小组就开发了GsTL，并将其开源，它仅包含头文件，使用起来很方便。
2、GsTL的总体设计
2.1 GsTL的组成
按照GSLIB的传统，GsTL同样也包含构建模块库和一系列程序模拟用于执行基本的地质统计学算法。但重点是构建模块库，使其能为未来软件开发提供整体框架。
- 基础库：这部分主要提供构建模块，使开发者能建立新程序实现新算法。地质统计估计和模拟常用的实体可以抽象为几个概念，每个概念可以定义为一系列类型，这些类型满足一些特定的要求。在c++中，这些概念通过类（class）表示，基础库也提供一些常用的地质统计学算法，Remy（2001）的文章有关于这些类和从用地质统计学算法的描述。
- 一系列程序：这部分的程序主要是空间变化性表征（变差函数模型、训练图像中结构提取等）、克里金估计、随机模拟。终端用户可以直接使用这些程序，不需任何修改。但是，这些程序模块必须重新组装，以满足他们的特定需求。
2.2 设计步骤
通过分析地质统计学算法和典型的应用问题寻找常用的概念，地质统计学算法常见的操作实体通过一系列概念表示。比如，很多地质统计估计和算法可以抽象为：
（1）建立网格系统，它包含一组位置，每个位置的属性值需要被估计或模拟；
（2）定义一个路径，可以访问到网格系统中的所有位置；
（3）对沿着路径的所有位置：
根据相邻数据评价该位置处的累计概率分布（ccdf）
从ccdf重采样赋值到所在位置
（4）结束
对这种通用的算法，我们可以定义一些重要的概念，包括：位置（location），属性值（property value），路径（path），网格系统（grid）,邻域（neighbourbood），累计概率分布古迹器（ccdf estimator），累计概率分布采样器（ccdf sampler）。这些概念可以是实体对象，比如位置，也可以是行为或函数，比如ccdf估计器。Remy2001对这些问题有详细讨论。
需要注意的是这些概念不是一小段c++代码，而是模块化成的c++类。这些类为维数相关信息提供了重要支撑，是对一组数据和算法的封装。
3 一些GsTL基础类型的设计与实现
对每个类型，首先展示声明函数，然后讨论它的设计和实现，紧接着是类型的模板参数，公开基类，它模拟的概念，它用到的类型及接口等。
3.1 点Point
template <class coordinate_type, size_t Dimension = 3>
struct Point
点是N维欧拉空间的元素，表达为笛卡尔坐标（p0, p1, p2, …, PN-1）,可以通过下标[]的方式访问点在某个维度的坐标，比如通过p[0]获取点p的第一个维度的坐标值。
除了支持+-*/操作，点的操作符+=, -=, -, ==, 已经重载，可以使用。<<和>>也重载了，方便输入和输出。
下面是个简单的实例
Point<int, 3> p = {{1, 2, 3}};
Point<int, 3> shift = {{1, 1, 1}};
Point<int, 3> q = p + shift; // q’s coordinates are (2, 3, 4)
q[0] = 5; // change q[0] from 2 to 5
3.2 地质值GeoValue
template<class location_type, class property_type>
class GeoValue
地质值是GeoValue位置和属性的点对，属性可以是连续变量，比如孔隙度、渗透率，也可以是离散变量，比如岩性。GeoValue也要一个标签指示是否为硬数据，已估计和未估计数据。GeoValue是很多地质统计学最基本的元素，因此它的设计对要率要求很高。
操作符<<和>>也支持，方便输出和输入
3.3 笛卡尔结构CartesianMesh
template<class location_type>
class CartesianMesh
笛卡尔网格通过初始点，每个维度上的维数，以及每个维度上的网格的大小组成。为了计算或搜索，很多地质统计学算法中地质值在网格中的存储都遵循一定的顺序。比如，在GISLIB中，地质值在规则的三维笛卡尔坐标系中通过网格索引排序：先一个一个点向东排序，然后一行一行向被排序，然后一层一层向上排序。非结构化的样本数据按照一维方式排序，通过索引访问。两种方式下都是通过笛卡尔网格构建地质网格的排序。因此，笛卡尔网格的职责就是翻译地质值的坐标与每个维度上的网格索引。
3.4 网格Grid
网格Grid是地质值GeoValue的集合，地质值不需要是均匀间隔，每个地质的位置都需要指定。网格内置了个笛卡尔网格坐标系，可以在构造函数中显示地输入，也可以通过网格中存储的地质值组成。
在搜索相邻网格过程中，需要在网格Grid中通过随机的方式访问地质值GeoValue，所以网格Grid也设计成了随机访问容器，即Grid提供了随机访问迭代器。
有时，我们需要迭代地访问网格中地质值的内容，比如搜索相邻数据计算两个位置的协方差，因此，也提供了带帮助函数的迭代器。
3.5 笛卡尔网格CartesianGrid
template <class geoValue_type>
class CartesianGrid : public Grid<geoValue_type>
可以看出笛卡尔网格继承自网格，它是均匀间隔的网格坐标系。在每个维度上通过索引即可计算坐标值。构建起来比Grid更简单。
3.6 网格适应器GridAdapter
template<class grid_type>
class GridAdapter
当利用笛卡尔网格的迭代器访问时，它先沿第一维度访问，然后第二维度，以此循环。这种方法对寻找地质值，搜索相邻数据很有优势。但是，我们优势需要按照不同的方式访问地质值，比如比如，很多序贯模拟算法要求按照随机的方式访问地质以避免人为性。
网格适应器接口可以用来提供一种获取相邻数据的方法，接口上洞的个数就是相邻数据的个数。从第一个观察洞，可以获取第一个最近的相邻数据，从第二个观察洞，可以获取第二临近的数据，依次类推。总之，网格适应器时细化的随机访问器，返回的地质值的索引。
3.7 搜索面积SearchArea
template<class location_type>
class SearchArea
很多地质统计学算法要求相邻网格的搜索要限制在一定范围内的子网格，这个区域成为搜索面域，它包含一系列网格的集合，这些网格分布在中心点的n维空间。它的返回值是围绕中心点的一组点。
3.8 邻域Neighborhood
template<class grid_type, class BinOp>
class Neighborhood : public GridAdapter<grid_type>
邻域的关键作用是为地质值寻找相邻数据提供服务，相邻数据的搜索要基于所复数的网格，并且搜索要限定在给定的面域内，用户需要指定计算两个地质值距离的方法。
邻域的搜索策略根据网格类型的不同分为两种：如果采用了笛卡尔网格，则采用螺旋式搜索方法；如果是非结构网格，则采用超级块搜索。
3.9 随机路径RandonPath
template<class grid_type>
class RandomPath : public GridAdapter<grid_type>
很多地质统计学算法都要求按照随机方法访问地质值，随机路径为地质值所在网格提供了接口，因此很容易实现序贯地访问网格节点。
因为网格适应器继承了随机路径，网格适应器已经实现了很多功能。除了构造器，随机路径为新的模拟路径提供了一种方法，它重新为访问节点提供随机路径，同时把网格标记为未模拟，除了已知的赢数据。
Reference
Zhanjun Ying. Design and implementation of some GsTL foundation classes. Stanford Center for Reservoir Forecasting, Department of Petroleum Engineering, Stanford University, CA 94305, April 25, 2001
Remy N , Shtuka A , Levy B , et al. GSTL: the geostatistical template library in C++[J]. Computers & Geosciences, 2002, 28(8):971-979.
