显示标签为“科学计算”的博文。显示所有博文
显示标签为“科学计算”的博文。显示所有博文

2009年2月28日星期六

GNU octave (Part I)

最近用了一次 octave 感觉很不错,和 Matlab 越来越像了,图形也做的很顺手了。加上它小巧(相对于 Matlab 而言),对 Matlab 的了解使我能很容易上手。顺便找了一下,manual 写的很详尽(我看的 3.1 版本的接近 600 pages),常用的功能都有。扫了一下主页,感觉还是很活跃的,于是下定决心好好研究一下,正好和 Matlab 比较促进我对 Matlab 的进一步了解。

就我观察,octave 不同版本号之间差异比较大,我现在说的 3.1 应该是和 Matlab 最接近的,而 2.x 还离得很远,这里我先就一些基本的问题说明一下。debian 里面 octave 一般分成下面几个 package,
  • octave*.*,这是 octave 的主程序(还有一个汇报 bug 的程序)所在,常用的函数一般有二进制版本(称为 built-in 函数)放在 /usr/lib/octave/version/oct/arch/*.oct。
  • octave*.*-common 里面主要是常用工具箱,m 文件多在 /usr/share/octave/version/m/toolbox-name/*.m。
  • octave*.*-doc 里面是文档,我现在就使用的 3.1 系列的文档。
  • octave*.*-htmldoc 还是文档。
下面这些类似于 toolbox 的东西来自 octave-forge 项目,
  • octave-ad 自动前向微商。
  • octave-audio 处理音频。
  • octave-bioinfo 处理生物信息学。
  • octave-combinatorics 组合学。
  • octave-comminiations 信号处理、通信。
  • octave-control 控制理论。
  • octave-data-smoothing 数据平滑化(似乎是用 Tikhonov regularizer?)。
  • octave-econometrics 度量经济学。
  • octave-epstk 生成 eps 文件。
  • octave-financial 金融工具箱。
  • octave-fixed 计算 fixed point。
  • octave-ftp 在 octave 里面使用 ftp。
  • octave-ga 就是遗传算法。
  • octave-general 广义函数。
  • octave-graceplot 原来 octave 调用 gnuplot 绘图,这个包用 grace。
  • octave-gsl 提供 GSL 的支持。
  • octave-ident 提供 system identification。
  • octave-image 提供图像处理的功能。
  • octave-informationtheory 提供信息论方面的函数。
  • octave-io 提供输入输出支持。
  • octave-irsa 提供 irregular sampling analysis。
  • octave-linear-algebra 线性代数。
  • octave-miscellaneous 一些杂七杂八的函数,如输出到 latex 等。
  • octave-missing-functions 搜索 octave 里面没有的 matlab 函数。
  • octave-nan 提供处理 NaN 的能力。
  • octave-nnet 前馈神经网络。
  • octave-netcdf 处理 NetCDF 数据。
  • octave-octgpr 散乱数据插值(看起来像 Gaussian process regression)。
  • octave-odebvp 线性常微分方程边值问题。
  • octave-odepkg 常微分方程初值问题。
  • octave-optim 优化工具箱。
  • octave-optiminterp 最优插值。
  • octave-outliers 发现 outlier。
  • octave-parallel 并行计算。
  • octave-pfstools 使用 octave 处理 HDR 图片。
  • octave-phisicalconstants 物理常量。
  • octave-plot 额外的作图工具。
  • octave-plplot 对 PLplot 库的支持。
  • octave-signal 信号处理。
  • octave-sockets 提供 socket 的接口。
  • octave-sp 半定规划。
  • octave-specfun 特殊函数。
  • octave-splines 样条。
  • octave-statistics 统计函数。
  • octave-strings 字符串处理。
  • octave-struct 对结构体的支持。
  • octave-symbolic 符号计算,基于 GiNaC 和 CLN。
  • octave-time 日期时间处理。
  • octave-tsa 时间序列分析。
  • octave-vrml 虚拟现实工具箱。
  • octave-senity 可能和 guide 类似。

以上可见 octave 并不比我们常用的 matlab 缺少很多功能,而且我们相信随着贡献者越来越多,一般科学研究需要的计算可以转交给 octave 来实现。我现在一些小的计算,懒的打开 Matlab,就会用 octave 来处理。数值计算的教学,我觉得有理由离开 Matlab 这个商业软件了。如果说 Scilab 并不顺手,那么 octave 已经不存在这个理由了。

不过,我们在这里还是列举一些 octave 和 Matlab 不同之处(参考这里),
  • 搜索路径应使用 addpath,而不是老版本直接对 LOADPATH 变量赋值。
  • Matlab 不支持 logical array 的 prod,但是 octave 支持。
  • Matlab 不支持对 nargin 赋值,但是 octave 支持。
  • Matlab 启动的时候执行 startup.m,octave 执行 /etc/octave*.conf。
  • Matlab v6 的 mat 文件可以被 octave 读取。v7 的可以被 2.9 之后的读取。
  • 不同长度的字符串在老版本的 octave 里面允许 vercat,而 matlab 不允许。
  • Matlab 允许 ! 开头表示系统调用,octave 只允许使用 system 函数。octave 里面 ! 和 ~ 都表示取反。
  • Matlab 不允许一个字符串和一个字符比较,但 octave 允许。
  • hist.m 函数在 octave 里面有一个 normalization,matlab 没有。
  • 老版本的 octave 通过环境变量 OCTAVE_VERSION 获得版本信息,matlab 和新版的 octave 使用 ver。
  • cell array 和 struct 的显示方是不同。
  • 函数 datastr 曾非常不一致。
  • Matlab 的 load 和 octave 的 load --force 等价。
  • Matlab 可以 load 读取空文件,octave 不允许。
  • Matlab 没有 printf 只有 fprintf 和 sprintf,octave 都有。
  • Matlab 的转置运算符和操作矩阵之间不允许有空格,octave 允许。
  • Matlab 一行分割成几行需要用 ...,octave 可以省略或者用 \ 也行。
  • 逻辑运算里面 Matlab 只允许 & 和 |,octave 允许 && ||(返回标量),和 & | (返回矩阵)。
  • Matlab 用 ^ 表示指数运算,octave 还允许用 **。
  • Matlab 的字符串用 ',octave 还可以用 "。
  • Matlab 结构控制语句都用 end 结束,octave 还允许 end if/for 等。
  • Matlab 调试 dbstep in 在 octave 为 dbstep、dbstep 在 octave 为 dbnext。
  • Matlab 广义特征值用 eig(A, B),在 octave 为 qz( A, B)。
  • Matlab 注释使用 %,octave 还可以使用 #。
  • Matlab 不允许对立即矩阵取 index,即 [ 1 2 3 4]( 3 ) 是不合法的,octave 允许。
  • Matlab 没有 += 这种类型的运算,octave 支持,并还支持 ++ --,和 C 一致。
  • Matlab 在新版本里面实现了 exception handling 和 octave 的 unwind protect 不同。
  • Matlab 实现 OOP 的风格不被 octave 接受。
  • Matlab 使用 private 子目录不被 octave 接受。

后面,我们开始更进一步学习 octave。

2008年7月30日星期三

GSL

GSL(GNU Scientific Library)是一个 C 写成的用于科学计算的库,下面是一些相关的包
Desired=Unknown/Install/Remove/Purge/Hold
| Status=Not/Inst/Cfg-files/Unpacked/Failed-cfg/Half-inst/trig-aWait/Trig-pend
|/ Err?=(none)/Hold/Reinst-required/X=both-problems (Status,Err: uppercase=bad)
||/ Name Version Description
+++-====================================-==========================-============================================
un gsl (no description available)
ii gsl-bin 1.11-2 GNU Scientific Library (GSL) -- binary packa
ii gsl-doc-pdf 1.11-2 GNU Scientific Library (GSL) Reference Manua
un gsl-ref-html (no description available)
un gsl-ref-pdf (no description available)
pn gsl-ref-psdoc (no description available)
un libgsl0 (no description available)
ii libgsl0ldbl 1.11-2 GNU Scientific Library (GSL) -- library pack
其中,libgsl0ldbl 是真正的库:
/usr/lib/libgsl.so.0.12.0
/usr/lib/libgslcblas.so.0.0.0
/usr/share/doc/libgsl0ldbl/changelog.Debian.gz
/usr/share/doc/libgsl0ldbl/changelog.gz
/usr/share/doc/libgsl0ldbl/THANKS.gz
/usr/share/doc/libgsl0ldbl/BUGS.gz
/usr/share/doc/libgsl0ldbl/README
/usr/share/doc/libgsl0ldbl/AUTHORS
/usr/share/doc/libgsl0ldbl/SUPPORT
/usr/share/doc/libgsl0ldbl/TODO.gz
/usr/share/doc/libgsl0ldbl/copyright
/usr/share/doc/libgsl0ldbl/NEWS.gz
/usr/share/lintian/overrides/libgsl0ldbl
/usr/lib/libgsl.so.0
/usr/lib/libgslcblas.so.0
另外还有一个相关的 gsl-bin,里面主要两个程序,
/usr/share/man/man1/gsl-histogram.1.gz
/usr/share/man/man1/gsl-randist.1.gz
/usr/share/doc
/usr/share/doc/gsl-bin
/usr/share/doc/gsl-bin/changelog.Debian.gz
/usr/share/doc/gsl-bin/changelog.gz
/usr/share/doc/gsl-bin/copyright
/usr/bin
/usr/bin/gsl-randist
/usr/bin/gsl-histogram
gsl-randist 和 gsl-histogram 分别是产生随机样本和计算直方图的程序。

我们这里根据 gsl-doc-pdf 给出如何使用 GSL 的程序接口和例程。GSL 的程序使用的头文件一般放在 /usr/include/gsl/ 目录(libgsl0-dev),C 程序通常用 #include 让预处理程序 cpp 读入对应的函数、宏声明,连接的时候通过 -lgsl -lgslcblas 到对应的库,通常还可能连接到 -lm。为了使用某些 inline 函数,可以打开 HAVE_INLINE 宏。注意由于 long double 和使用平台相关,一般不推荐使用。

GSL 里面的命名规则大致是使用 gsl 作为前缀(没有 namespace,sigh),函数的话一般是 gsl_foo_fn(对应 double)其余的为 gsl_foo_type_fn,对于类型一般是 gsl_foo 或者 gsl_foo_type(没有 template,sigh again)。对应的头文件一般是 gsl_foo.h(含有所有的类型)或者 gsl_foo_type.h。

GSL 里面出错处理遵循 POSIX 线程库,正常返回 0,出错返回非零,并依照 gsl_errno.h 里面设置出错值。可以用 gsl_strerror 将返回值用字符串表达出来。默认情况下 GSL 提供的 error handler 就是打印出错问题并调用 abort,这是一个 gsl_error_handler_t 类型的函数,可以通过 gsl_set_error_handler() 函数设定。

Mathematical Functions 常用数学函数 gsl_math.h
包括常用的数学常数(M_*),GSL_POSINF、GSL_NEGINF、GSL_NAN 以及相应判断的函数 gsl_isnan()、gsl_isinf() 和 gsl_finite()。另外提供了一些函数值快速计算的方法,gsl_log1p() 计算 log( 1 + x),gsl_exp1m() 计算 ex-1,gsl_hypot() 和 gsl_hypot3 计算欧氏空间范数,gsl_acosh()、
gsl_asinh()、gsl_atanh() 是反双曲函数,gsl_ldexp(x, y) 计算 x . 2y,gsl_frexp() 计算在二进制下科学记数法下 x 的基数部分。求幂次 gsl_pow_int 或者 gsl_pow_n(n = 2, ..., 9)。测试符号 GSL_SIGN,奇偶性 GSL_IS_EVEN 和 GSL_IS_ODD。取大小 GSL_MAX、GSL_MIN。浮点数大小最好用 gsl_fcmp 函数。

Complex Numbers 复数 gsl_complex.h、gsl_complex_math.h
定义一个复数可以用 gsl_complex_rect 或者 gsl_complex_polar,另外获得实部、虚部 GSL_REAL 和 GSL_IMAG,设定 GSL_SET_COMPLEX,GSL_SET_REAL 和 GSL_SET_IMAG。可以求辐角 gsl_complex_arg()、模长 gsl_complex_abs()、模长平方 gsl_complex_abs2、模长对数 gsl_complex_logabs()。复数的加减乘除就是 gsl_complex_op() 其中 op 可以为 add、sub、mul 和 div,另外和实数有类似的运算 gsl_complex_op_real(),和虚数有 gsl_complex_op_imag(),共轭 gsl_complex_conjugate()、逆 gsl_complex_inverse() 和相反数 gsl_complex_negative。另外如求平方根 gsl_complex_sqrt(对实数求加 _real),幂次 gsl_complex_pow(次数为实数加 _real),指数 gsl_complex_exp,对数 gsl_complex_log(loh10 或者 log_b)。另外还有三角函数、反三角函数、双曲函数、反双曲函数。

Polynomial 多项式 gsl_poly.h
一般可以给定多项式的系数,用一个数组从低阶到高阶,调用 Horner 法求多项式函数值可以用 gsl_poly_eval(),对复数而言是 gsl_poly_complex_eval(),对复系数多项式为 gsl_complex_poly_complex_eval()。另外一种表达方式是使用 Newton 差分法表达,这时输入的是差分节点创建一个多项式,gsl_poly_dd_init(),求函数值可用 gsl_poly_dd_eval(),还可以把这种类型的多项式转换成为 Taylor 展开的形式 gsl_poly_dd_taylor()。对二次多项式求解根可以用 gsl_poly_solve_quadratic(),复根可以用 gsl_poly_complex_solve_quadratic()。对三次方程用 cubic 替换 quadratic。对高于 4 次的多项式没有解析解,往往用矩阵特征值进行近似,GSL 提供了一种解法,首先用 gsl_poly_complex_workspace_alloc() 分配储存根的空间,然后调用 gsl_poly_complex_solve() 求解,求解之后应该用 gsl_poly_complex_workspace_free() 释放。

Special Functions 特殊函数 gsl_sf.h、gsl_sf_*.h
一般有两种调用形式,一种和正常函数类似直接 gsl_sf_function(),另一种是 gsl_sf_function_e() 将返回值的地址传入函数。gsl_sf_result.h 提供了一个用于估计误差的结构体,一般函数有三个 mode 控制计算精度 GSL_PREC_DOUBLE、GSL_PREC_SINGLE 和 GSL_PREC_APPROX。提供的特殊函数有 airy(见 gsl_sf_airy.h)的函数值、零点、导数、导数零点,Bessel 函数(见 gsl_sf_bessel.h)的函数值、零点,Clausen 函数(见 gsl_sf_clausen.h),Coulomb 函数(见 gsl_sf_coulomb.h),Coupling 系数(见 gsl_sf_coupling.h),Dawson 函数(见 gsl_sf_dawson.h),Debye 函数(见 gsl_sf_debye.h),Dilogorithm 函数(见 gsl_sf_dilog.h),乘法误差函数(见 gsl_sf_elementary.h),椭圆积分(见 gsl_sf_ellint.h),Jacobi 椭圆函数(见 gsl_sf_elljac.h),误差函数(见 gsl_sf_erf.h,GNU libc 也有类似的函数),指数函数(见 gsl_sf_exp.h),指数积分(见 gsl_sf_expint.h),Fermi-Dirac 函数(见 gsl_sf_fermi_dirac.h),Gamma 和 Beta 函数(见 gsl_sf_gamma.h),Gegenbauer 函数(见 gsl_sf_gengenbauer.h),超几何函数(见 gsl_sf_hyperg.h),Laguerre 函数(见 gsl_sf_laguerre.h),Lambert W 函数(见 gsl_sf_lambert.h),Legendre 函数和球面调和函数(见 gsl_sf_legendre.h),对数及其相关函数(见 gsl_sf_log.h),Mathieu 函数(见 gsl_sf_mathieu.h),幂函数(见 gsl_sf_pow_int.h),Psi(digamma) 函数(见 gsl_sf_psi.h),Synchrotron 函数(见 gsl_sf_synchrotron.h),transport 函数(见 gsl_sf_transport.h),三角双曲函数(见 gsl_sf_trig.h),Zeta 函数(见 gsl_sf_zeta.h)。

Vectors and Matrices 向量与矩阵 gsl_block.h、gsl_vector.h、gsl_matrix.h
创建 vector 也好、matrix 也好,都依赖 gsl_block 这个结构,可以用 gsl_block_alloc() 和 gsl_block_calloc() 分配,gsl_block_free() 释放,另外有对流的输入输出,如 gsl_block_fread()、gsl_block_fwrite()、gsl_block_fprintf() 和 gsl_block_fscanf()。不管是 vector 还是 matrix 都含有一个 gsl_block 的指针,操作和 block 类似。如 vector 的的类似操作就是 gsl_vector_alloc()、
gsl_vector_calloc()、gsl_vector_fread()、gsl_vector_fwrite()、gsl_vector_fprintf() 和 gsl_vector_fscanf(),另外可以通过 gsl_vector_get() 和 gsl_vector_set() 获得/设定某一分量的值,gsl_vector_ptr() 和 gsl_vector_const_ptr() 获得一分量的地址,另外有一些函数方便初始 vector,如 gsl_vector_set_all()、gsl_vector_set_zero() 和 gsl_vector_set_basis()。为了访问一个 vector 元素的子集,可以使用 vector view 对象,这可以用一些函数产生(有对应的 const 版本),并最好仅仅在 stack 内使用(也就是直接操作对象本身,而不是指针),如 gsl_vector_subvector() 产生一个连续的子集,gsl_vector_subvector_with_stride() 产生一个带固定间隔的子集,gsl_vector_complex_real() 和 gsl_vector_complex_imag() 产生一个 real 或者 image 部分的 view,gsl_vector_view_array() 对一个数组产生 vector view,gsl_vector_view_array_with_stride() 产生带有固定间隔的 vector view。vector 之间的复制或者互换有 gsl_vector_memcpy() 和 gsl_vector_swap()。vector 元素之间互换 gsl_vector_swap_elements(),逆序 gsl_vector_reverse()。vector 之间的四则运算 gsl_vector_op(),数乘(op = scale),加上常数(op=add_constant)。vector 最大最小(op=max、min、minmax)或者对应的 index(op = max_index、min_index、minmax_index)。判断一个 vector 是否为 0 向量(op=isnull),正(ispos),负(isneg),非负(isnonneg)。matrix 和 vector 稍微不同之处在于用两个下标索引,前面函数多数只要把 vector 换成 matrix 即可。另外还可以为 matrix 的行或者列建立 view,gsl_matrix_(sub)row/column(),或者对角元素 gsl_matrix_(sub, super)diagonal()。将矩阵一行/列读到/写到一个 vector 可以用 gsl_matrix_get/set_row/col()。矩阵行列互换 gsl_matrix_swap_rows/columns() 或者方阵的行列交换 gsl_matrix_swap_rowcol(),转置或转置复制 gsl_matrix_transpose()、gsl_matrix_transpose_memcpy()。矩阵运算中 mul_elements 和 div_elements 是对元素的。

Permutations 置换 gsl_permutation.h
这是产生置换的基本数据结构,一般用 gsl_permutation_(c)alloc() 分配内存,gsl_permutation_init() 初始化为置换幺元,可以用 gsl_permutation_memcpy() 进行复制,gsl_permutation_free() 释放。访问置换元素可以用 gsl_permutation_get(),互换用 gsl_permutation_swap()。另外可以用 gsl_permutation_size() 获得置换大小,gsl_permutation_data() 获得指向置换的指针,gsl_permutation_valid() 验证是否为合法的 permutation。另外有一些置换的操作,如 gsl_permutation_reverse() 逆转,gsl_permutation_inverse() 求逆,依照字典序计算下一个/前一个有 gsl_permutation_next/prev()。将 permutation 应用到数组上可以用 gsl_permute,或者逆 psl_permute_inverse(),对 vector 可以 gsl_permute_vector(_inverse)(),几个置换可以相乘 gsl_permutation_mul()。类似的 permutation 也有输入输出函数。另外置换存在一种正则表达方式可以用 gsl_permutation_linear_to_canonical() 转换,可以计算一个 permutation 含有几个 cycle 等。

Combinations 组合 gsl_combination.h
和置换类似的结构,但是处理的是组合问题。

Sorting 排序 gsl_heapsort.h、gsl_sort_*.h
首先提供了一个 quick sort 的补充的 heapsort,gsl_heapsort() 和 gsl_heapsort_index()。排序数组 or vector 可以用 gsl_sort() 或者 gsl_sort_vector(),另外也有带索引的版本。求最小/大的 k 个元素,可以用 gsl_sort(_vector)_smallest/largest(_index)()。

BLAS Support 基本线性代数子程序支持 gsl_blas.h、gsl_cblas.h
BLAS 支持三个 level 的运算,level 1 是 vector 的,level 2 是 matrix-vector 的,level 3 是 matrix-matrix 的操作。操作对象的类型为 SDCZ 对应 float、double、float complex 和 double complex,矩阵的特性为 GE(一般)、GB(一般带状矩阵)、SY(对称)、SB(对称带状)、SP(对称,packed)、HE、HB、HP(Hermite)、TR、TB、TP(三角阵)。操作类型有 DOT(内积)、AXPY(a x + y)、MV(矩阵 x 向量)、SV(矩阵逆乘向量)、MM(矩阵相乘)、SM(矩阵的逆乘另外一个矩阵)。GSL 提供的命令形式为 gsl_blas_*。

Linear Algebra 线性代数 gsl_linalg.h
这部分包括了最常用的数值线性代数运算,如矩阵 LU 分解求解线性方程组(带 permutation,可以 inplace 等版本),QR 分解(包括选取列的),SVD,Cholesky 分解,实对称矩阵的对角化分解(本征分解),Hermite 矩阵的对角化分解,实矩阵的 Hessenberg 分解,实矩阵对的 Hessenberg 分解,双对角化(bidiagonalization),Householder 变换,Householder 变换求解线性方程组,三对角阵,balancing(通过相似变换使得行列的范数相当)。

Eigensystems 求解特征值 gsl_eigen.h
这部分包括实对称矩阵的特征值、Hermite 矩阵的特征值,以及非对称矩阵的特征值(利用 Schur 分解)以及对应的广义特征值问题求解的函数。一般需要 alloc 一个 workspace,然后调用对应的函数计算特征值、特征向量,最后 free 掉 workspace。另外还提供了对应的函数用于同时整理特征值与特征向量。

Fast Fourier Transform 快速 Fourier 变换 gsl_fft_*.h
快速 Fourier 变换这里分成了对复数、实数(更困难一些,需要保证逆变换获得是实数,使用的也是 halfcomplex 表达系数)两种处理。而对于数据为 2 的幂次的可以直接用 Cooley-Tuckey 算法,不是的话另外有一套算法(需要预先分配 workspace)。每套算法提供 forward(计算 Fourier 变换),inverse(逆变换),backward(不带规范化常数的逆变换)和 transform(通过参数选择 forward 还是 backward)。

Numerical Integration 数值积分 gsl_integration.h
函数命名方式是 gsl_integration_*(),Q 表示 quadrature routine,N 和 A(表示是否自适应),G 和 W(一般积分和带权值函数的积分),S 和 P(容易消解的奇点或者提供特别困难的点),I(无穷积分),O(振荡积分),F(Fourier 积分),C(Cauchy 主值)。积分里面设置停止条件是设置相对误差或者绝对误差。

Random Number Generation 随机数生成器 gsl_rsg.h
首先需要 gsl_rng_alloc() 一个对应的类型,然后 gsl_rng_set() 设置 seed,gsl_rng_free() 释放。同时还可以通过环境变量 GSL_RNG_TYPE 和 GSL_RNG_SEED 以及函数 gsl_rng_env_setup() 获取,然后通过设定对应的生成器就可以利用 gsl_rng_uniform() 产生 [0, 1) 的均匀分布,gsl_rng_uniform_pos() 产生 (0, 1) 的均匀分布,以及 gsl_rng_uniform_int() 产生指定范围内的均匀整数分布。另外可以通过 gsl_rng_name() 获得该生成器的名称,gsl_rng_get() 返回一个在 gsl_rng_min() 和
gsl_rng_max() 之间的随机数。如果需要更细致的处理生成器,还提供了一些函数用于处理它的状态 IO。另有一章详细介绍各种分布下随机数的生成情况。

Quasi Random Sequences 拟随机序列 gsl_qrng.h
与前一章不同的是不需要初始 seed,调用结构和前一章类似。

Random Number Distribution 随机数分布 gsl_randist.h
这里包含了绝大多数常用分布,命名规则如下 gsl_ran_dist(_),这里 dist 是分布名称,如 gaussian 等,后如果没有 _ 表示产生随机数,如果 _pdf 是密度,分布函数用了两种形式,一般是 cdf_dist_P 和 cdf_dist_Q 以及对应的逆 Pinv 和 Qinv。一共有下面几种分布,gaussian、gaussian_tail、bivariate_gaussian、exponential、laplace、exppow、cauchy、rayleigh、rayleigh_tail、landau、levy、levy_skew、gamma、flat、lognormal、chisq、fdist、tdist、beta、logistic、pareto、dir_2d、weibull、gumbel1、gumbel2、dirichlet。对于离散分布,有限个取值可以用 gsl_ran_discrete_preproc() 将分布密度列(或者差一个 scale factor)转换成为一个 gsl_ran_discrete_t 类型的结构,并传递给 gsl_ran_discrete() 产生随机数,gsl_ran_discrete_pdf() 产生分布列,产生的结构可以用 gsl_ran_discrete_free() 释放。另外提供了 poisson、bernoulli、binomial、multinomial、negative_binomial、pascal、geometric、hypergeometric、logarithmic。除了分布函数,还可以把指定序列随机打乱 gsl_ran_shuffle()、随机选元素 gsl_ran_choose()、选择一个子集 gsl_ran_sample()。

Statistics 统计 gsl_stats.h
主要提供统计函数,如求均值 gsl_stats_mean(),子样方差(无偏) gsl_stats_variance()、子样方差(已知期望,有偏) gsl_stats_variance_m(),标准差两个版本 std 和 std_m,与期望平方和 tss 和 tss_m,另外两个是 variance_with_fixed_mean 和 sd_with_fixed_mean。绝对偏差 absdev 和 absdev_m,skew 用 skew、skew_m_sd,峰度 kurtosis 和 kurtosis_m_sd,自相关性 lag1_autocorrelation 和 lag1_autocorrelation_m,协方差 covariance 和 covariance_m,相关系数 correlation,另外有对应的带权值版本在前面加上 w 即可。另外也提供了最大最小以及对应 index 的函数,计算中位数以及分位数的函数。

Histograms 直方图 gsl_hostogram*.h
分一维和两维,差异不大,主要把一维的增加新的一个变量,命名在 histogram 后面加 2d 即可。大致使用的方法是,首先 gsl_histogram_alloc() 分配空间,然后 gsl_histogram_set_ranges 设置节点(或者用 gsl_histogram_set_ranges_uniform() 设置均匀的节点),最后 gsl_histogram_free()。另外还提供了复制 gsl_histogram_memcpy() 以及自身复制 gsl_histogram_clone() 的函数。可以通过 gsl_histogram_increment() 增加元素到计数,也可以 gsl_histogram_accumulate() 增加任意权值(计数器用的实数),可以获得某个 bin 的权值 gsl_histogram_get(),或者某个 bin 上下界 gsl_histogram_get_range(),整个 histogram 的上下界 gsl_histogram_min/max(),bin 的数目 gsl_histogram_bins(),而 gsl_histogram_reset() 将整个 histogram 清零。gsl_histogram_find() 返回某个值所在的 bin。另有 max/min_val/bin 返回最大值或者出现的 bin,还有利用这个 histogram 计算 mean、sigma(标准差)、sum。另外两个 histogram 可以用 gsl_histogram_equal_bins_p() 看看是否可以使用 add/sub/mul/div 等运算,shift 可以让所有值 + 常数,scale 是乘。还有一些 IO 的函数。可以用 gsl_histogram 来创建一个 gsl_histogram_pdf,这个就和前面讲的随机数一样用。

N-tuples N 元组 gsl_ntuple.h
非常简单的数据结构,用于把数据写入/读出文件,提供的基本操作有 gsl_ntuple_create() 创建一个空文件(截断已存在),gsl_ntuple_open() 打开已存在文件,或者 gsl_ntuple_write() 将 ntuple 写入文件,或者从文件中读入 ntuple 数据 gsl_ntuple_read(),最后需要 gsl_ntuple_close() 关闭。可以把一个 ntuple 数据读入喂给 histogram 进行统计,这主要使用 gsl_ntuple_project() 函数。

Monte Carlo Integration 蒙特卡罗积分 gsl_monte_*.h
实现的是最基本的三种积分方法,在 gsl_monte.h 里面声明了积分函数的基本形式 gsl_monte_function,在 gsl_monte_plain.h 提供的是最简单均匀采样积分方法,首先 gsl_monte_plain_alloc() 分配 workspace,gsl_monte_plain_init() 初始化,然后 gsl_monte_plain_integrate() 积分,最后 gsl_monte_plain_free() 释放 workspace。在 gsl_monte_miser.h 里面使用的分层蒙特卡罗积分,可以用 gsl_monte_miser_state 结构控制算法细节。而 gsl_monte_vegas.h 是使用 impartance sampling,gsl_monte_vegas_state 控制算法细节。复杂的 MCMC 等这里并没有实现。

Simulated Annealing 模拟退火 gsl_siman.h
只有一个函数 gsl_siman_solve(),提供优化函数等信息即可使用。

Ordinary Differential Equations 常微分方程 gsl_odeiv.h
主要使用 gsl_odeiv_system 结构,需要提供方程标准形式的函数和偏导(即 Jacobi 矩阵)。另外和算法相关的是 stepping function,如 Runge-Kutta 方法等,也有自适应版本的,这种函数的目的是为了计算指定一个 step 下函数值。GSL 还提供了计算一个区间内函数变化(若干 step)的函数(evolve)。但是由于对常微分方程数值解不很了解,这里就略去吧。

Interpolation 插值 gsl_spline.h、gsl_interp.h
提供了三次样条和 Akima 样条。比较 low-level 的函数为用户提供了非常细致的控制,通过 gsl_interp_alloc() 分配需要的空间,选择适当的算法,gsl_interp_init() 初始化节点,最后 gsl_interp_free() 释放。为了搜索某个位置(以便计算函数值)可以用 gsl_interp_bsearch(),也可以使用 gsl_interp_accel 对象(先 init,然后 find 和 free)。另外提供了最常用的函数值、一阶导、二阶导和积分以及对应使用 gsl_interp_accel 的接口。high-level 的函数主要在 gsl_spline.h 里面提供。和 gsl_interp_* 系列相似。

Numerical Difference 差分 gsl_deriv.h
提供了中心差分 gsl_derive_central()、前向 gsl_deriv_forward()、后向 gsl_deriv_backward()。

Chebyshev Approximation 车比雪夫逼近 gsl_chebyshev.h
提供了 [-1, 1] 上一组正交多项式,这对应的是 1/sqrt(1-x^2) 为权的函数空间。首先用 gsl_cheb_alloc() 分配空间产生 gsl_cheb_series 结构,然后 gsl_cheb_init(),最后 gsl_cheb_free()。提供了计算函数值、函数值误差,导函数和积分。

Series Acceleration 级数加速 gsl_sum.h
提供了一个 Levin u-transform 的东西,没听说过。其意义在于减少求和项,提高计算精度。使用方式就是通过 gsl_sum_levin_u_alloc() 分配 workspace,继而通过 gsl_sum_levin_u_accel() 分配,最后 gsl_sum_levin_free()。如果不需要估计误差,则可以更快。

Wavelet Transform 小波变换 gsl_wavelet*.h
与 FFT 类似,但是没有 backward 类型,分一维和二维,有 daubechies、haar 和 bspline。

Discrete Hankel Transform 离散汉克尔变换 gsl_dht.h
与 FFT 类似,但是是对极坐标的,调用方式和 FFT 类似。

One-dimensional Root Finding 一维函数求零点 gsl_roots.h
有两种方式(基于搜索 gsl_root_fsolver 和基于导数 gsl_root_fdfsolver)。首先选取合适的 solver,命名方式都是 gsl_root_*solver_type,然后 alloc。之后可以用 gsl_root_*solver_set() 设定初始态,开始迭代使用 gsl_root_*solver_iterate(),也可以直接用 gsl_root_*solver_root() 求出根。另可以 gsl_root_fsolver_x_upper() 和 gsl_root_fsolver_x_lower() 返回控制根的区间。通过 gsl_root_test_*() 可以测试相对误差、残差。提供的算法有 bisection、falsepos 和 brent,利用梯度的有 newton、secant、steffenson。

One-dimensional Minimization 一维函数求极小 gsl_min.h
最小化在某种意义上就是求导函数的零点。因此调用方法和前一章极为类似。算法有 goldensection 和 brent。

Multi-dimensional Root Finding 多维函数求零点 gsl_multiroots.h
类似。算法有 hybridsj,hybridj,newton,gnewton。不用梯度的算法有 hybrids,hybrid,dnewton,broyden。

Multi-dimensional Minimization 多维函数求极小 gsl_multimin.h
类似。算法有 conjugate_fr、conjugate_pr、bfgs、bfgs2、steepest_descent、nmsimplex。

Least Square Fitting 最小二乘拟合 gsl_fit.h
分单变量和多变量。gsl_fit_linear() 和 gsl_fit_wlinear() 分别是线性和加权线性问题的单变量拟合(即线型回归),另外 gsl_fit_linear_est() 还估计误差。多元情形为 gsl_fit_mul()、gsl_fit_wmul() 与 gsl_fit_mul_linear()。对广义的 LSF 问题,需要使用 gsl_multifit_linear_alloc() 分配 workspace,最后释放,类似的函数有 gsl_multifit_linear 和 weighted 版本,另有 _svd 版本,使用 SVD 计算结果。

Nonlinear Least Square Fitting 非线性最小二乘拟合 gsl_multifit_nlin.h
与多维函数最小化类似。

Basic Splines 基础样条 gsl_bspline.h
先 gsl_bspline_alloc() 产生 workspace,然后 gsl_bspline_knots() 或者 gsl_bspline_knots_uniform() 设置节点,gsl_bspline_eval() 计算函数值,最后 gsl_bspline_free() 释放空间。

Physical Constants 物理常数 gsl_const_mksa.h
各种可能用到的物理常数,命名一般为 GSL_CONST_MKSA_*。

IEEE Floating Point Arithmetic 浮点算术 gsl_ieee_utils.h
提供了输出 float 和 double 的 gsl_ieee_printf_*,另外可以用 gsl_ieee_env_setup() 设置对应的计算环境。

相关程序:
  • GSL 的 C++ wrapper:GSLwrap
  • 一个用 GSL 写的 PCA 程序,和 mex 一起工作。