GNU Scientific Library examples rewritten

  • Difficulty: 020/100
  • Author: jeremy theler
  • Keywords: VAR, VECTOR, VECTORS, SIZE, DATA, NUMBER, PRINT, SEP, PRINT_FUNCTION, FORMAT, TEXT, SEP, SEPARATOR, PARAMETRIC, MIN, MAX, TYPE, OUTER_STEPS, PHASE_SPACE, INTERPOLATION, MINIMIZE, ALGORITHM, GUESS, GRADIENT, GRADTOL, TOL, IF, ENDIF, FIT, TO, VIA, COLUMNS, FILE_PATH, VERBOSE, integral, random, .=, interpolation, derivative, root, sqrt, func_min, cos, random, sobol, vecmax, vecmin, cspline, cspline_periodic, conjugate_fr, nmsimplex2, done_outer, step_inner, step_outer, in_static, integral,

As wasora heavily relies on the “GNU Scientific Library”: “http://www.gnu.org/software/gsl/” to perform most of its numerical-related tasks, it should be possible to run most of GSL’s examples with a wasora input. Indeed, this section shows how some of the GSL examples found in the documentation can be implemented by just writing a wasora input file instead of coding an compiling and ad-hoc C program. However, as wasora is not (and does not want to be) a programming language, some of these examples may look awkward but it should be kept in mind that they are written as to resemble as much as possible the original GSL examples written in C.

The motley examples are taken from the “GNU Scientific Library v1.16 documentation”: “http://www.gnu.org/software/gsl/manual/html_node/”, and are sorted in increasing complexity (and not by its appearance order in the manual). Each case refers to the original section of the manual where it is discussed. First, the C code from the official GSL distribution is shown. Then, the corresponding wasora input that should accomplish the same computation is listed and finally a terminal where the problem is solved in both ways. If the case involve the generation of figures, these are included afterwards.

The paragraphs below are a combination of texts that belong to the original GSL documentation with explanations written by the author of wasora. Comments and extensions are welcome.

integration.was

Consider the following integral, which has an algebraic-logarithmic singularity at the original

!bt [ _0^1 x^{-1/2} (x) , dx = 4 ] !et

The original example uses the QAGS integrator, which is not implemented in wasora. However, the regular integral functional (which corresponds to the QAG integrator) does an acceptable job. Currently, wasora does not provide a way to report the estimated errors.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/integration.c

The corresponding wasora input:

# GSL example 17.14, numerical integration
VAR x
PRINT %.18f integral(x^(-1/2)*log(x),x,0,1)
$ gcc -o integration integration.c -lgslcblas -lgsl -lm
$ ./integration
result          = -4.000000000000085265
exact result    = -4.000000000000000000
estimated error =  0.000000000000135447
actual error    = -0.000000000000085265
intervals =  8
$ wasora integration.was
-3.999889192698361295
$ 

diff.was

This example estimates the derivative of the function \(f(x) = x^{3/2}\) at \(x=2\) and at \(x=0\). The function \(f(x)\) is undefined for \(x<0\) so the derivative at \(x=0\) is computed using gsl_deriv_forward in the C example. To reproduce this behavior in wasora, a positive value \(p>0\) is passed as the optional fifth argument (see the derivative functional documentation in the wasora reference). Again, wasora does not provide a way to estimate the error.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/diff.c

The corresponding wasora input:

# GSL example 28.2, differentiation

f(x) := x^(3/2)

PRINT TEXT "f(x) = x^(3/2)"
PRINT TEXT "x = 2.0"
PRINT TEXT "f'(x) = " %.10f derivative(f(x),x,2,1e-8) SEP " "
PRINT TEXT "exact = " %.10f 1.5*sqrt(2) SEP " "
PRINT
PRINT TEXT "x = 0.0"
PRINT TEXT "f'(x) = " %.10f derivative(f(x),x,0,1e-8,1) SEP " "
PRINT TEXT "exact = " %.10f 0 SEP " "
$ gcc -o diff diff.c -lgslcblas -lgsl -lm
$ ./diff
f(x) = x^(3/2)
x = 2.0
f'(x) = 2.1213203120 +/- 0.0000005006
exact = 2.1213203436

x = 0.0
f'(x) = 0.0000000160 +/- 0.0000000339
exact = 0.0000000000
$ wasora diff.was
f(x) = x^(3/2)  
x = 2.0 
f'(x) =  2.1213203120
exact =  2.1213203436

x = 0.0 
f'(x) =  0.0000000160
exact =  0.0000000000
$ 

You may want to test what happens if the derivative at \(x=0\) is computed without the optional fifth argument \(p\).

roots.was

This example uses the function solver gsl_root_fsolver_brent for Brent’s method to solve the following equation

!bt [ x^2 - 5 = 0 ] !et

with solution \(x = \sqrt{5} = 2.236068\dots\). First, the function to be solved is defined in the following two files demo_fn.h and demo_fn.c:

@@@CODE 020-gsl/demo_fn.h

@@@CODE 020-gsl/demo_fn.c

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/roots.c

The corresponding wasora input:

# GSL example 33.10, root finding
a = 1
b = 0
c = -5
f(x) := (a*x + b)*x + c 

PRINT %.7f root(f(x),x,0,5) sqrt(5)
$ gcc -o roots roots.c -lgslcblas -lgsl -lm
$ ./roots
using brent method
 iter [    lower,     upper]      root        err  err(est)
    1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000
    2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000
    3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000
    4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000
    5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300
Converged:
    6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666
$ wasora roots.was
2.2360680   2.2360680
$ 

As an illustration, the wasora version (which is significantly smaller than its corresponding C partner) uses the default tolerance for the root finding routines instead of the \(0.001\) value of the original example, hence the difference in the results.

min.was

The following program uses the Brent algorithm to find the minimum of the function \(f(x) = \cos(x) + 1\), which occurs at \(x = \pi\). The starting interval is \([0,6]\). The GSL version uses an initial guess for the minimum of \(x=2\). However, wasora takes the initial guess at the center of the starting interval.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/min.c

The corresponding wasora input:

# GSL example 34.8, one-dimensional minimization
VAR x
PRINT %.7f func_min(cos(x)+1,x,0,6)
$ gcc -o min min.c -lgslcblas -lgsl -lm
$ ./min
using brent method
 iter [    lower,     upper]       min        err  err(est)
    0 [0.0000000, 6.0000000] 2.0000000 -1.1415927 6.0000000
    1 [2.0000000, 6.0000000] 3.5278640 +0.3862713 4.0000000
    2 [2.0000000, 3.5278640] 3.1748217 +0.0332290 1.5278640
    3 [2.0000000, 3.1748217] 3.1264576 -0.0151351 1.1748217
    4 [3.1264576, 3.1748217] 3.1414743 -0.0001183 0.0483641
    5 [3.1414743, 3.1748217] 3.1415930 +0.0000004 0.0333474
Converged:
    6 [3.1414743, 3.1415930] 3.1415927 +0.0000000 0.0001187
$ wasora min.was
3.1415914
$ 

As seen in the terminal, the functional func_min performs the minimization at once and does not provide a mean to show information about intermediate steps.

rngunif.was

The following program demonstrates the use of a random number generator to produce uniform random numbers in the range \([0.0, 1.0)\). In the GSL version, The numbers depend on the seed used by the generator. The default seed can be changed with the GSL_RNG_SEED environment variable to produce a different stream of numbers. The generator itself can be changed using the environment variable GSL_RNG_TYPE. In wasora, the generator is fixed to the gsl_rng_knuthran2002 method, basically because the author looks up at Knuth as the most important figure of modern computer science. Therefore, to obtain two identical sequences of numbers both with GSL and wasora, we run both examples with the same generator and seed.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/rngunif.c

The corresponding wasora input:

# GSL example 18.13, random number generation
static_steps = 10
PRINT %.5f random(0,1,$1)
$ gcc -o rngunif rngunif.c -lgslcblas -lgsl -lm
$ GSL_RNG_SEED=123 GSL_RNG_TYPE=knuthran2002 ./rngunif
GSL_RNG_TYPE=knuthran2002
GSL_RNG_SEED=123
0.32585
0.53345
0.27282
0.84846
0.64330
0.90798
0.19579
0.63717
0.92833
0.71629
$ wasora rngunif.was 123
0.32585
0.53345
0.27282
0.84846
0.64330
0.90798
0.19579
0.63717
0.92833
0.71629
$ 

qrng.was

Quasi-random sequences are used by wasora to perform parametric studies. The idea is to sweep a certain portion of the parameter space by choosing test points, which can be linear (type LINEAR), logarithmic (LOGARITHMIC), uniformly chosen at random (RANDOM), randomly distributed using a gaussian distribution (GAUSSIANRANDOM) or follow a certain quasi-random sequence (types SOBOL, NIEDERREITER, HALTON and REVERSEHALTON). By defining a two-dimensional window \([0,1]\times[0,1]\) and just printing the chosen parameters, the GSL example can be reproduced in wasora.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/stat.c

The corresponding wasora input:

# GSL example 19.6 quasi-random sequences
PARAMETRIC x y MIN 0 0 MAX 1 1 OUTER_STEPS 1024 TYPE sobol
PRINT %.5f x y SEP " "
$ gcc -o qrng qrng.c -lgslcblas -lgsl -lm
$ ./qrng > qrng-gsl.dat
$ wasora qrng.was | tee qrng-was.dat | graph -Tpng -C -m -3 -S 2 > qrng.png
$ diff -s qrng-gsl.dat qrng-was.dat
Files qrng-gsl.dat and qrng-was.dat are identical
$ 
qrng.png
qrng.png

stat.was

A basic example of how to use the statistical functions.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/stat.c

The corresponding wasora input:

# GSL example 21.10, statistics
VECTOR data SIZE 5 DATA 17.2 18.1 16.5 18.3 12.6

N = vecsize(data)
mean = vecsum(data)/N
variance = sum((data(i) - mean)^2/(N-1), i, 1, N)

PRINT TEXT "the dataset is" %g data SEPARATOR ", "
PRINT TEXT "the sample mean is       " %g mean
PRINT TEXT "the estimated variance is" %g variance
PRINT TEXT "the largest value is     " %g vecmax(data)
PRINT TEXT "the smallest value is    " %g vecmin(data)
$ gcc -o stat stat.c -lgslcblas -lgsl -lm
$ ./stat
The dataset is 17.2, 18.1, 16.5, 18.3, 12.6
The sample mean is 16.54
The estimated variance is 5.373
The largest value is 18.3
The smallest value is 12.6
$ wasora stat.was
the dataset is, 17.2, 18.1, 16.5, 18.3, 12.6
the sample mean is          16.54
the estimated variance is   5.373
the largest value is        18.3
the smallest value is       12.6
$ 

ode-initval.was

The following program solves the second-order nonlinear Van der Pol oscillator equation

!bt [ u’‘(t) + u’(t) (u(t)^2 - 1) + u(t) = 0 ] !et

This can be converted into a first order system suitable for use with the routines described in this chapter by introducing a separate variable for the velocity, \(v = u'(t)\),

!bt \begin{align*} u' &= v\cr v' &= -u + \mu v (1-u^2) \end{align*}

!et

The program begins by defining functions for these derivatives and their Jacobian. The main function uses driver level functions to solve the problem. The program evolves the solution from \((u, v) = (1,0)\) at \(t=0\) to \(t=100\). The step-size \(h\) is automatically adjusted by the controller to maintain an absolute accuracy of \(10^{-6}\) in the function values \((u, v)\).
The loop in the example prints the solution at the points \(t_i = 1, 2, \dots, 100\).

It is important to note that wasora does not use GSL’s ordinary differential equations solvers because wasora works with DAEs instead of ODEs (which of course are a particular case of DAEs). The solver is based on the SUNDIALS library IDA. Nevertheless, dynamical systems are one of the main subjects of wasora so this example is particularly relevant. The special variable rel_error is set to \(10^{-6}\) as in GSL and both the variables min_dt and max_dt are set to one, so as to reproduce the solution only at the points \(t_i=1,2,\dots,100\).

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/ode-initval.c

The corresponding wasora input:

# GSL example 26.6, ordinary differential equations
# (solved with sundials, not with gsl)
PHASE_SPACE u v

# solver settings
rel_error = 1e-6
end_time = 100
min_dt = 1   # report results at t=0,1,2,...,100
max_dt = 1

# dynamical system parameter
mu = 10

# initial conditions
u_0 = 1
v_0 = 0

# dynamical system equations
u_dot .= v
v_dot .= -u + mu*v*(1-u^2) 

PRINT t u v
$ gcc -o ode-initval ode-initval.c -lgslcblas -lgsl -lm
$ ./ode-initval > vanderpoolgsl.dat
$ wasora ode-initval.was > vanderpoolwas.dat
$ pyxplot vanderpool.ppl
$ 
vanderpool.png
vanderpool.png

interp.was

The following program demonstrates the use of the interpolation and spline functions. It computes a cubic spline interpolation of the 10-point dataset \((x_i, y_i)\) where \(x_i = i + \sin(i)/2\) and \(y_i = i + \cos(i^2)\) for \(i = 0 \dots 9\).

To keep the original spirit, the data is plotted using GNU ploutils graph progam.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/interp.c

The corresponding wasora input:

# GSL example 27.7.1, interpolation
n = 10
VECTOR fx SIZE n
VECTOR fy SIZE n
FUNCTION yi(x) VECTORS fx fy INTERPOLATION cspline

# remember the example has C indexes!
fx(i) = (i-1) + sin(i-1)/2
fy(i) = (i-1) + cos((i-1)^2)

# write the original data
PRINT TEXT "\#m=-3,S=3"
PRINT_FUNCTION yi FORMAT %g
# print the interpolated data
PRINT TEXT "\#m=3,S=0"
PRINT_FUNCTION yi MIN fx(1) MAX fx(n) STEP (fx(n)-fx(1))/100 FORMAT %g
$ gcc -o interp interp.c -lgslcblas -lgsl -lm
$ ./interp > interpgsl.dat
$ head interpgsl.dat
#m=0,S=2
0 1
1.42074 1.5403
2.45465 1.34636
3.07056 2.08887
3.6216 3.04234
4.52054 5.9912
5.86029 5.87204
7.32849 7.30059
8.49468 8.39186
$ tail interpgsl.dat
9.11 9.56698
9.12 9.58874
9.13 9.61052
9.14 9.63232
9.15 9.65414
9.16 9.67598
9.17 9.69783
9.18 9.71969
9.19 9.74156
9.2 9.76343
$ graph -Tpng -C -L GSL < interpgsl.dat > interpgsl.png
$ wasora interp.was > interpwas.dat
$ head interpwas.dat
#m=-3,S=3   
0   1
1.42074 1.5403
2.45465 1.34636
3.07056 2.08887
3.6216  3.04234
4.52054 5.9912
5.86029 5.87204
7.32849 7.30059
8.49468 8.39186
$ tail interpwas.dat
8.37751 8.23336
8.46957 8.35571
8.56164 8.49437
8.6537  8.6486
8.74576 8.8162
8.83782 8.99493
8.92988 9.18258
9.02194 9.3769
9.114   9.57568
9.20606 9.77669
$ graph -Tpng -C -L wasora < interpwas.dat > interpwas.png
$ 
interpgsl.png
interpgsl.png
interpwas.png
interpwas.png

interpp.was

The next program demonstrates a periodic cubic spline with 4 data points. Note that the first and last points must be supplied with the same \(y\)-value for a periodic spline.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/interpp.c

The corresponding wasora input:

# GSL example 27.7.2, interpolation
FUNCTION f(x) INTERPOLATION cspline_periodic DATA {
0.00  0.15
0.10  0.70
0.27 -0.10
0.30  0.15
}

# original data
PRINT TEXT "\#m=-3,S=3"
PRINT_FUNCTION f FORMAT %g
# interpolated data
PRINT TEXT "\#m=3,S=0"
PRINT_FUNCTION f MIN f_a MAX f_b STEP (f_b-f_a)/100 FORMAT %g
$ gcc -o interpp interpp.c -lgslcblas -lgsl -lm
$ ./interpp > interppgsl.dat
$ graph -Tpng -C -L GSL < interppgsl.dat > interppgsl.png
$ wasora interpp.was > interppwas.dat
$ graph -Tpng -C -L wasora < interppwas.dat > interppwas.png
$ 
interppgsl.png
interppgsl.png
interppwas.png
interppwas.png

multimin.was

This example program finds the minimum of a paraboloid function. The location of the minimum is offset from the origin in \(x\) and \(y\), and the function value at the minimum is non-zero. The initial step-size is chosen as 0.01, a conservative estimate in this case, and the line minimization parameter is set at 0.0001. The program terminates when the norm of the gradient has been reduced below 0.001.

In wasora, as the minimization function is algebraic, the gradient can be entered easily as shown. If it was not entered, it would be automatically computed by finite differences (new versions of GSL do the same).

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/multimin.c

The corresponding wasora input:

# GSL example 36.9.1, multidimensional minimization
a = 10
b = 20
c = 30
x0 = 1
y0 = 2
f(x,y) := a*(x-x0)^2 + b*(y-y0)^2 + c

MINIMIZE f ALGORITHM conjugate_fr {
     GUESS 5          7
  GRADIENT 2*a*(x-x0) 2*b*(y-y0)
      STEP 1e-2
   GRADTOL 1e-3
       TOL 1e-4   
}

IF done_outer
 PRINT  "\# minimum found at"
ENDIF
IF step_inner=1|done_outer
 PRINT  %5.0f step_outer %.5f x y %10.5f f(x,y) SEPARATOR " "
ENDIF
$ gcc -o multimin multimin.c -lgslcblas -lgsl -lm
$ ./multimin
    1 4.99629 6.99072  687.84780
    2 4.98886 6.97215  683.55456
    3 4.97400 6.93501  675.01278
    4 4.94429 6.86073  658.10798
    5 4.88487 6.71217  625.01340
    6 4.76602 6.41506  561.68440
    7 4.52833 5.82083  446.46694
    8 4.05295 4.63238  261.79422
    9 3.10219 2.25548   75.49762
   10 2.85185 1.62963   67.03704
   11 2.19088 1.76182   45.31640
   12 0.86892 2.02622   30.18555
Minimum found at:
   13 1.00000 2.00000   30.00000
$ wasora multimin.was
    0 0.00000 0.00000  120.00000
    1 4.99629 6.99072  687.84780
    2 4.98886 6.97215  683.55456
    3 4.97400 6.93501  675.01278
    4 4.94429 6.86073  658.10798
    5 4.88487 6.71217  625.01340
    6 4.76602 6.41506  561.68440
    7 4.52833 5.82083  446.46694
    8 4.05295 4.63238  261.79422
    9 3.10219 2.25548   75.49762
   10 1.20067 -2.49832  435.09973
   11 2.19088 1.76182   45.31640
   12 0.86892 2.02622   30.18555
   13 3.51283 1.49743   98.19448
# minimum found at  
   13 1.00000 2.00000   30.00000
$ 

This example shows an important difference. The GSL program prints a line only after successful steps of the optimization method. But wasora prints a line with the instantaneous position each time the objective function \(f(x,y)\) has to be computed, as the whole input is executed to evaluate the function to be minimized (this example is simple, but wasora was designed to optimize complex computations such as nuclear reactor designs). Therefore, GSL shadows the real number of evaluations (which are directly related to the computational cost of the problem). Analyze with care.

nmsimplex.was

Here is another example using the Nelder-Mead Simplex algorithm to minimize the same example object function, as above.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/nmsimplex.c

The corresponding wasora input:

# GSL example 36.9.2, multidimensional minimization
a = 10
b = 20
c = 30
x0 = 1
y0 = 2
f(x,y) := a*(x-x0)^2 + b*(y-y0)^2 + c

MINIMIZE f ALGORITHM nmsimplex2 {
     GUESS 5          7
  GRADIENT 2*a*(x-x0) 2*b*(y-y0)
      STEP 1          1
       TOL 1e-2
}

IF in_outer_initial
 PRINT  TEXT "\#    x        y         z"
ENDIF
IF done_outer
 PRINT  TEXT "\# converged to minimum at"
ENDIF
IF (not(in_outer_initial)&step_inner=1)|done_outer
 PRINT  %5.0f step_outer %10.3e x y TEXT "f() =" %7.3f f(x,y) SEPARATOR " "
ENDIF
$ gcc -o nmsimplex nmsimplex.c -lgslcblas -lgsl -lm
$ ./nmsimplex
    1  6.500e+00  5.000e+00 f() = 512.500 size = 1.130
    2  5.250e+00  4.000e+00 f() = 290.625 size = 1.409
    3  5.250e+00  4.000e+00 f() = 290.625 size = 1.409
    4  5.500e+00  1.000e+00 f() = 252.500 size = 1.409
    5  2.625e+00  3.500e+00 f() = 101.406 size = 1.847
    6  2.625e+00  3.500e+00 f() = 101.406 size = 1.847
    7 -1.776e-15  3.000e+00 f() =  60.000 size = 1.847
    8  2.094e+00  1.875e+00 f() =  42.275 size = 1.321
    9  2.578e-01  1.906e+00 f() =  35.684 size = 1.069
   10  5.879e-01  2.445e+00 f() =  35.664 size = 0.841
   11  1.258e+00  2.025e+00 f() =  30.680 size = 0.476
   12  1.258e+00  2.025e+00 f() =  30.680 size = 0.367
   13  1.093e+00  1.849e+00 f() =  30.539 size = 0.300
   14  8.830e-01  2.004e+00 f() =  30.137 size = 0.172
   15  8.830e-01  2.004e+00 f() =  30.137 size = 0.126
   16  9.582e-01  2.060e+00 f() =  30.090 size = 0.106
   17  1.022e+00  2.004e+00 f() =  30.005 size = 0.063
   18  1.022e+00  2.004e+00 f() =  30.005 size = 0.043
   19  1.022e+00  2.004e+00 f() =  30.005 size = 0.043
   20  1.022e+00  2.004e+00 f() =  30.005 size = 0.027
   21  1.022e+00  2.004e+00 f() =  30.005 size = 0.022
   22  9.920e-01  1.997e+00 f() =  30.001 size = 0.016
   23  9.920e-01  1.997e+00 f() =  30.001 size = 0.013
converged to minimum at
   24  9.920e-01  1.997e+00 f() =  30.001 size = 0.008
$ wasora nmsimplex.was
#    x        y         z   
    1  6.000e+00  6.000e+00 f() = 600.000
    2  5.500e+00  5.000e+00 f() = 412.500
    3  6.750e+00  2.000e+00 f() = 360.625
    4  5.500e+00  1.000e+00 f() = 252.500
    5  4.000e+00  3.000e+00 f() = 140.000
    6  2.875e+00  5.000e-01 f() = 110.156
    7 -1.776e-15  3.000e+00 f() =  60.000
    8 -2.500e-01  6.000e+00 f() = 365.625
    9 -5.313e-01  1.375e+00 f() =  61.260
   10  2.352e+00  7.812e-01 f() =  77.974
   11 -1.248e+00  2.477e+00 f() =  85.079
   12  1.588e+00  2.564e+00 f() =  39.834
   13  1.261e+00  1.651e+00 f() =  33.118
   14  1.760e+00  1.804e+00 f() =  36.551
   15  7.173e-01  1.828e+00 f() =  31.390
   16  9.134e-01  2.131e+00 f() =  30.417
   17  7.181e-01  2.088e+00 f() =  30.951
   18  1.097e+00  2.060e+00 f() =  30.167
   19  1.000e+00  1.962e+00 f() =  30.029
   20  1.085e+00  1.948e+00 f() =  30.127
   21  9.954e-01  2.043e+00 f() =  30.037
   22  1.047e+00  1.986e+00 f() =  30.026
   23  1.015e+00  2.019e+00 f() =  30.009
   24  9.732e-01  1.984e+00 f() =  30.012
# converged to minimum at   
   24  9.920e-01  1.997e+00 f() =  30.001
$ 

nlfit.was

The following example program fits a weighted exponential model with background to experimental (simulated) data , \(Y = A \exp(-\lambda t) + b\).

The wasora implementation, instead of generating the random gaussian data (which can be nevertheless implemented with the random_gauss function), the data generated by the GSL program is sent to a file which is used to define the objective function. This situation on the one hand is closer to actual applications of wasora to fit real experimental data and, on the other hand, we make sure we use the same set of data in both cases so we expect to find the same results.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/nlfit.c

The corresponding wasora input:

# GSL example 38.9, non-linear fit
#
# GRADIENT is optional, wasora will compute it numerically if not provided
# VERBOSE activates step-by-step information
# 

# initial guess
A = 0
lambda = 1
b = 0
Y(t) := A * exp(-lambda*t) + b

FUNCTION data(t) FILE_PATH experimental.dat COLUMNS 2 3

FIT Y TO data VIA A lambda b {
  DELTAEPSREL 1e-4 DELTAEPSABS 1e-4
#   GRADIENT exp(-lambda*t) -t*A*exp(-lambda*t) 1
  VERBOSE
}

IF done_outer
 PRINT %.5f TEXT "A      =" A      TEXT "+/-" sigma_A      SEP " " 
 PRINT %.5f TEXT "lambda =" lambda TEXT "+/-" sigma_lambda SEP " "
 PRINT %.5f TEXT "b      =" b      TEXT "+/-" sigma_b      SEP " "
ENDIF

IF done_outer
 PRINT_FUNCTION  FILE_PATH nlfit.dat data Y
ENDIF
$ gcc -o nlfit nlfit.c -lgslcblas -lgsl -lm
$ ./nlfit | tee gsl.out
data: 0 6.01339 0.1
data: 1 5.51538 0.1
data: 2 5.26109 0.1
data: 3 4.77746 0.1
data: 4 4.45135 0.1
data: 5 3.9049 0.1
data: 6 3.50439 0.1
data: 7 3.415 0.1
data: 8 3.24274 0.1
data: 9 3.1222 0.1
data: 10 2.83763 0.1
data: 11 2.5347 0.1
data: 12 2.43917 0.1
data: 13 2.38083 0.1
data: 14 2.31609 0.1
data: 15 2.06083 0.1
data: 16 1.94568 0.1
data: 17 1.91413 0.1
data: 18 1.75951 0.1
data: 19 1.66507 0.1
data: 20 1.73793 0.1
data: 21 1.57552 0.1
data: 22 1.52507 0.1
data: 23 1.40961 0.1
data: 24 1.39521 0.1
data: 25 1.41689 0.1
data: 26 1.37604 0.1
data: 27 1.26095 0.1
data: 28 1.28963 0.1
data: 29 1.42267 0.1
data: 30 1.22829 0.1
data: 31 1.19918 0.1
data: 32 1.18999 0.1
data: 33 0.930076 0.1
data: 34 1.22461 0.1
data: 35 1.14738 0.1
data: 36 1.114 0.1
data: 37 1.19512 0.1
data: 38 1.26958 0.1
data: 39 1.06198 0.1
iter:   0 x =      1.00000000      0.00000000      0.00000000 |f(x)| = 117.349
status = success
iter:   1 x =      1.64656260      0.01814555      0.64656260 |f(x)| = 76.4594
status = success
iter:   2 x =      2.85864925      0.08091443      1.44791355 |f(x)| = 37.6859
status = success
iter:   3 x =      4.94892814      0.11943460      1.09463075 |f(x)| = 9.58212
status = success
iter:   4 x =      5.02174105      0.10287695      1.03389255 |f(x)| = 5.63076
status = success
iter:   5 x =      5.04520407      0.10405524      1.01941629 |f(x)| = 5.44398
status = success
iter:   6 x =      5.04535782      0.10404906      1.01924871 |f(x)| = 5.44397
chisq/dof = 0.800996
A      = 5.04536 +/- 0.06028
lambda = 0.10405 +/- 0.00316
b      = 1.01925 +/- 0.03782
status = success
$ cat gsl.out | grep data > experimental.dat
$ wasora nlfit.was
fititeration: 0  0.000e+00  1.000e+00  0.000e+00    status = success    |r|^2 = 16.7263
fititeration: 1  5.040e+00  1.000e+00  2.076e+00    status = success    |r|^2 = 6.6945
fititeration: 2  5.033e+00  8.766e-01  2.075e+00    status = success    |r|^2 = 6.49716
fititeration: 3  5.021e+00  6.777e-01  2.070e+00    status = success    |r|^2 = 6.05048
fititeration: 4  4.986e+00  4.101e-01  2.043e+00    status = success    |r|^2 = 4.95003
fititeration: 5  4.686e+00  1.104e-01  1.872e+00    status = success    |r|^2 = 4.53229
fititeration: 6  4.805e+00  1.449e-01  1.557e+00    status = success    |r|^2 = 1.80061
fititeration: 7  4.922e+00  9.566e-02  1.091e+00    status = success    |r|^2 = 1.06197
fititeration: 8  5.036e+00  1.041e-01  1.027e+00    status = success    |r|^2 = 0.545438
fititeration: 9  5.045e+00  1.040e-01  1.019e+00    status = success    |r|^2 = 0.544395
fititeration: 10     5.045e+00  1.040e-01  1.019e+00    status = success    |r|^2 = 0.544395
A      = 5.04536 +/- 0.60279
lambda = 0.10405 +/- 0.03157
b      = 1.01925 +/- 0.37824
# chisq/dof = 0.00800988
# A = 5.0454 # +/- 0.60279
# lambda    = 0.10405 # +/- 0.031571
# b = 1.0192 # +/- 0.37824
$ pyxplot nlfit.ppl
$ 
nlfit.png
nlfit.png

Note that wasora cannot currently handle uncertainties on the experimental data, so they are assumed to be equal to one (instead of the GSL example where they are identically equal to 0.1 for every point). Therefore, the reported uncertainties on the coefficients are larger in wasora than in the GSL example.