Optiminterp

From OCG Test Wiki

1 Optimal interpolation Fortran module with Octave interface

[edit] Optimal interpolation Fortran module with Octave interface

This Fortran 90 module performs a n-dimensional optimal interpolation (OI). The optimal interpolation allows to interpolate arbitrarily located observations to a regular grid using a background field as first guess. The observations and the background fields may contain errors. Optimal interpolation merges observations and background taking their expected variances into account. The merged field is optimal in the sense that it has the lowest error variance.

The error correlation function of the background field is assumed to decrease exponentially with the square of the distance. The correlation length in every dimension must be specified. Other correlation functions can easily be implemented.

In some disciplines, this technique is called objective analysis or kriging but the underling algorithm is the same.

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

This package is also included in the octave-forge repository. Some Linux distributions already include octave-forge and you only need to install the octave-forge package in order to use this package. For other distribution, you can compile either octave-forge or optiminterp from source.

[edit] Implementation

The vector $\mathbf x^b$ is the background state estimation, a first guess. The observation are contained in the vector $\mathbf y^o$ . The operator $\mathbf H$ extracts from a state vector the corresponding values at the location of the observation. The vector $\mathbf H \mathbf x^b$ represents thus the first guess at the location of the observations. The difference is noted $\mathbf d$ :

$\mathbf d =\mathbf y^o - \mathbf H \mathbf x^b$

The elements of matrix $\mathbf P^b$ the represents to covariance between two points $(x 1,..., x n)$ and $(y 1,..., y n)$ in the n-dimensional space. Knowing the error covariance of the observations $\mathbf R$ and the error covariance of the background state $\mathbf P^b$ the Kalman gain matrix is defined as:

$\mathbf K =\mathbf P^b \mathbf H^T \left( \mathbf H \mathbf P^b \mathbf H^T+ \mathbf R \right)^{-1}$

The optimal state vector $\mathbf x^a$ can be computed as:

$\mathbf x^a - \mathbf x^b = \mathbf K \mathbf d$

The optimal interpolation code operates on difference from the background state. From the difference $\mathbf d$ it computes the state vector correction $\mathbf x^a -\mathbf x^b$ . It is thus up to the users of this program to form the vector $\mathbf d$ and after state vector correction is obtained to compute the analysis $\mathbf x^a$ . \\

The error covariance $\mathbf P^a$ of the analyzed state $\mathbf x^a$ is given by:

$\mathbf P^a = \mathbf P^b - \mathbf K \mathbf H \mathbf P^b$

The following assumptions are made:

$\mathbf R$ is assumed diagonal.
$\mathbf P^b$ is simply given by the following parametrization:

$P^b(x_1,...,x_n,y_1,...,y_n)=\sigma(x_1,...,x_n)^2 \exp\left( -\frac{(x_1-y_1)^2}{L_1^2} ... -\frac{(x_n-y_n)^2}{L_n^2} \right)$

The error variance $σ(x 1,..., x n) 2$ and the correlation lengths $L 1$ ,..., $L n$ have to be specified.

Only a limited number $m$ (specified to the user) contribute to the analysis of a given point.

[edit] Requirements

Fortran 90 compiler
LAPACK and BLAS

For the optional GNU Octave interface you will need:

Octave with its development package

[edit] Download the package

Download the sources of the toolbox from http://ocgmod1.marine.usf.edu/OI/optiminterp-0.2.1.tar.gz.

wget http://ocgmod1.marine.usf.edu/OI/optiminterp-0.2.1.tar.gz

[edit] Fortran interface

Decompress the file

tar -xzvf optiminterp-0.2.1.tar.gz
cd optiminterp-0.2.1

The optimal interpolation module is written in Fortran 90. The file example_optiminterp.F90 in src shows how to use the optimal interpolation module from Fortran. On a Linux or UNIX system, you would compile this example, by

cd src
gfortran -o example_optiminterp optimal_interpolation.F90  example_optiminterp.F90 -llapack -lblas

Replace gfortran with the name of your Fortran 90 compiler. You can check, if the example is correctly compiled by:

$ ./example_optiminterp
Expected:  0.022206
Obtained:  0.022206

[edit] Octave interface

[edit] Automatic installation using Octave's package manager

If you have octave 2.9.9 or newer and if octave was compiled with a Fortran 90 compiler, you can use Octave's package manager to compile and install the optimal interpolation toolbox. In the directory which contains the file optiminterp-0.2.1.tar.gz, start a octave session and enter the following commands, to install and load and package:

octave:1> pkg install optiminterp-0.2.1.tar.gz
octave:2> pkg load all

You should see no error messages. If you are curious to know where octave has placed the optimal interpolation toolbox issue the following command:

octave:3> pkg list
Package Name | Version | Installation directory
-------------+---------+-----------------------
 optiminterp |   0.2.1 | /home/abarth/octave//optiminterp-0.2.1

If your version of octave was compiled with g77 or with a problematic gfortran (see troubleshooting below), you can tell mkoctfile to use a different Fortran compiler. To use the g95 compiler, you have to set the following environment variables before starting octave and installing the package.

export F77=g95
export FFLAGS=-fno-second-underscore
export FLIBS=$(g95 -print-file-name=libf95.a)

[edit] Manual installation

Decompress the file:

tar -xzvf optiminterp-0.2.1.tar.gz
cd optiminterp-0.2.1/src

Since we need to link Fortran and C++ code, the manual installation is compiler specific. Basically we need to know where the Fortran run-time libraries are located. The locations and the names of the run-time libraries are specified with the FLIBS environment variable.

[edit] gfortran

If you have gfortran version 4.1 or newer, all you need to do is to compile with:

gfortran -fPIC -c optimal_interpolation.F90 optiminterp_wrapper.F90
FLIBS=$(gfortran -print-file-name=libgfortran.so) mkoctfile optiminterp.cc optiminterp_wrapper.o optimal_interpolation.o

Unfortunately, the Fortran 90 interpolation module doesn't work with gfortran 4.0. It compiles but it generates a run-time error. In this case, I recommend to use the g95 Fortran compiler.

[edit] g95

As gfortran, g95 is also an open source Fortran 90 compiler. Binaries for most platforms are available at the g95 website. Download and extract the tar file corresponding to your platform, e.g.:

cd $HOME
wget -O - http://ftp.g95.org/g95-x86_64-32-linux.tgz | tar  xvfz -

This installs g95 in your home directory under g95-install. In order to avoid some typing you can define an alias. You need probably to adapt the name of the g95 executable as it depends on your platform.

alias g95="$HOME/g95-install/bin/x86_64-unknown-linux-gnu-g95"

This alias is only defined in the current shell. In order to make this alias permanent you need to add this line in your $HOME/.bashrc file.

Issue the following command:

g95  -fno-second-underscore  -fPIC -c optimal_interpolation.F90 optiminterp_wrapper.F90
FLIBS=$(g95 -print-file-name=libf95.a) mkoctfile optiminterp.cc optiminterp_wrapper.o optimal_interpolation.o

[edit] Intel compiler

Assuming that your Intel compiler is installed in /opt/intel/fc/9.0/, you can compile it with:

ifort -fPIC -O3 -c optimal_interpolation.F90 optiminterp_wrapper.F90
FLIBS="-L/opt/intel/fc/9.0/lib -lifcore_pic  -lguide" mkoctfile -v optiminterp.cc optiminterp_wrapper.o optimal_interpolation.o

[edit] PGI compiler

The PGI compiler works also:

pgf90 -fPIC -c optimal_interpolation.F90 optiminterp_wrapper.F90
DL_LD="pgf90 -fPIC" mkoctfile -v optiminterp.cc optiminterp_wrapper.o optimal_interpolation.o

or assuming that the Fortran run-time libraries are in /usr/pgi/linux86-64/6.0/libso/, you can compile it with:

pgf90 -fPIC -c optimal_interpolation.F90 optiminterp_wrapper.F90
FLIBS="-L/usr/pgi/linux86-64/6.0/libso/ -lpgf90 -lpgf90_rpm1 -lpgf902 -lpgf90rtl -lpgftnrtl -lpgc" mkoctfile -v optiminterp.cc optiminterp_wrapper.o optimal_interpolation.o

[edit] Testing the optimal interpolation module

The octave script "test_optiminterp.m" performs a 1D, 2D, 3D and 4D optimal interpolation. All tests should pass; any error indicates that either there is a bug in the optimal interpolation package or that it is incorrectly installed.

octave:1> test_optiminterp
Testing 1D-optimal interpolation: OK
Testing 2D-optimal interpolation: OK
Testing 3D-optimal interpolation: OK
Testing 4D-optimal interpolation: OK

[edit] Using the optimal interpolation module

The following is a minimal example about how to use the optimal interpolation package in Octave. This example is included in the file "example_optiminterp.m". The package also includes an example about how to use the OI package from Fortran 90 called "example_optiminterp.f".

% number of observations to interpolate

on = 200;

% create randomly located observations within 
% the square [0 1] x [0 1] 

x = rand(1,on);
y = rand(1,on); 

% the underlying function to interpolate

f = sin(6*x) .* cos(6*y);

% the error variance of the observations

var = 0.01 * ones(on,1);

% the grid onto which the observations are interpolated

[xi,yi] = ndgrid(linspace(0,1,100));

% the correlation length in x and y direction

lenx = 0.1;
leny = 0.1;

% number of influential observations

m = 30;

% run the optimal interpolation
% fi is the interpolated field and vari is its error variance

[fi,vari] = optiminterp2(x,y,f,var,lenx,leny,m,xi,yi);

[edit] Compilation with OpenMP

The optimal interpolation module is parallelized using OpenMP directives. To use multiple CPUs of a SMP machine you need to enable OpenMP (using the -openmp flag of the Intel compiler or -mp flag of PGI compiler) and set the environment variable OMP_NUM_THREADS to the desired number of CPUs to use before starting octave.

[edit] Troubleshooting

Undefined symbols with Intel's Fortran Compiler.

octave:1> test_optiminterp
Testing 1D-optimal interpolation: failed
error: [...]/optiminterp-0.2.1/optiminterp.oct: undefined symbol: __intel_cpu_indicator

You will also need to link against Intel's ifport library by adding "-lifport" at the end of the link command. If you get other error messages complaining about undefined symbols you might want to follow the instructions described in Intel's Compilers for Linux - Compatibility with the GNU Compilers (Thanks Claudio for pointing this out)

Runtime error with gfortran.

octave:1> test_optiminterp
Testing 1D-optimal interpolation: Fortran runtime error: dimension of return array incorrect

This error was observed with "gcc version 4.1.2 20060613 (prerelease) (Debian 4.1.1-5)" and with "gcc (GCC) 4.0.2 20051125 (Red Hat 4.0.2-8)" on Fedora Core 4. This is probably a bug in gfortran. This problem does not occur in more recent versions of gfortran. The latest gfortran binaries are available at http://quatramaran.ens.fr/~coudert/gfortran/.

[edit] Comments and feedback

Any suggestions, comments, bug fixes, ... are of course very appreciated. Send them to: abarth at marine dot usf dot edu.

Alexander Barth, last update January, 11 2007.

Retrieved from "http://ocgmod1.marine.usf.edu/wiki/index.php/Optiminterp"

Navigation