Tutorial

Installation

This package is supported just for Julia version 1.0. Consequently, it uses package 3.0. Currently RAFF is registered in General Julia Registers, so the package can be installed using the Julia package manager. From the Julia REPL, type ] to enter into Pkg REPL mode and run:

pkg> add RAFF

In what follows, we provide some simple examples on how to solve problems with RAFF. All the examples, and some other ones, are given in the examples/ directory as Julia scripts.

Basic usage

Just to illustrate the potential and basic usage of RAFF, let us consider the following data set given by an array $A$:

Let's suppose the first column of $A$ as an experimental measure with result given by second column. It is easy to see that the fitting function in this case is accurate and given by

Now let's perturb one result of second column of $A$. For example, consider $A_{6,2} = 2.55$. Assuming the model for fitting given by

we have by classical least squares as result θ = [0.904329, 1.3039]. On the other hand, when we consider RAFF algorithm we obtain the correct answer θ = [1.0, 1.0]. Moreover, we also have a list of the possible outliers.

In order to run RAFF algorithm we need to setup

julia> using RAFF

and define the data set and model:

julia> A=[-2.0  5.0;
          -1.5  3.25;
          -1.0  2.0 ;
          -0.5  1.25;
           0.0  1.0 ;
           0.5  2.55;
           1.0  2.0 ;
           1.5  3.25;
           2.0  5.0 ;];

julia> model(x, θ) = θ[1] * x[1]^2 + θ[2]
model (generic function with 1 method)

After that, we can run method raff:

julia> raff(model, A, 2)
** RAFFOutput **
Status (.status) = 1
Solution (.solution) = [0.999999, 1.0]
Number of iterations (.iter) = 39
Number of trust points (.p) = 8
Objective function value (.f) = 3.672268261041584e-11
Number of function evaluations (.nf) = 39
Number of Jacobian evaluations (.nj) = 39
Index of outliers (.outliers) = [6]

The number 2 above is the number of variables in model, i. e., the number of parameters to adjust in the model. The output is a RAFFOutput type. For example, to access only the parameters of solution, we can use

julia> output = raff(model, A, 2)
** RAFFOutput **
Status (.status) = 1
Solution (.solution) = [0.999999, 1.0]
Number of iterations (.iter) = 39
Number of trust points (.p) = 8
Objective function value (.f) = 3.672268261041584e-11
Number of function evaluations (.nf) = 39
Number of Jacobian evaluations (.nj) = 39
Index of outliers (.outliers) = [6]

julia> output.solution
2-element Array{Float64,1}:
 0.999998774479423
 1.000003442880683

Note that RAFF algorithm detects and ignores possible outliers. In order to see which points are outliers, we can access the outliers attribute.

julia> output.outliers
1-element Array{Int64,1}:
 6

More details about RAFFOutput type and other options can be obtained in API section.

By default RAFF uses automatic differentiation, more specifically ForwardDiff.jl package. But is possible to call RAFF methods with gradient vector of model. For example, considering the above example, we have,

Programming this gradient and running RAFF we have

julia> gmodel!(g, x, θ) = begin
          g[1] = x[1]^2
          g[2] = 1.0
       end
gmodel! (generic function with 1 method)

julia> raff(model, gmodel!, A, 2)
** RAFFOutput **
Status (.status) = 1
Solution (.solution) = [0.999999, 1.0]
Number of iterations (.iter) = 39
Number of trust points (.p) = 8
Objective function value (.f) = 3.672268261041584e-11
Number of function evaluations (.nf) = 39
Number of Jacobian evaluations (.nj) = 39
Index of outliers (.outliers) = [6]

Preliminary tests have shown that the use of explicit derivatives is 10 times faster than automatic differentiation.

Multivariate models

RAFF supports the use of multivariate fitting functions to data sets of different dimensions. To illustrate how this works, consider the following example:

julia> data = [1.0 1.0    2.0
               0.0 0.0    4.0
               7.0 1.5   -4.5
               2.0 2.0  -17.0 # outlier
               0.0 8.6   -4.6]
5×3 Array{Float64,2}:
 1.0  1.0    2.0
 0.0  0.0    4.0
 7.0  1.5   -4.5
 2.0  2.0  -17.0
 0.0  8.6   -4.6

and the following model

julia> model(x, θ) = θ[1] * x[1] + θ[2] * x[2] + θ[3]
model (generic function with 1 method)

Note that this model has two variables $(x_1, x_2)$ and three parameters $(\theta_1, \theta_2, \theta_3)$. This problem has one outlier (data[4,:]), so there are 4 trusted points. Let's run RAFF and check the answer.

julia> output = raff(model, data, 3)
** RAFFOutput **
Status (.status) = 1
Solution (.solution) = [-0.999999, -0.999999, 4.0]
Number of iterations (.iter) = 28
Number of trust points (.p) = 4
Objective function value (.f) = 2.9377653725356374e-11
Number of function evaluations (.nf) = 28
Number of Jacobian evaluations (.nj) = 28
Index of outliers (.outliers) = [4]

The right answer is [-1.0, -1.0, 4.0]. As we can note, RAFF get a good fit for the data set. Handling the output follows the same pattern as the one-dimensional case.

In order to get improvements in processing time, we can code the gradient vector of model too:

julia> gmodel!(g, x, θ) = begin
           g[1] = x[1]
           g[2] = x[2]
           g[3] = 1.0
       end
gmodel! (generic function with 1 method)

julia> output = raff(model, gmodel!, data, 3)
** RAFFOutput **
Status (.status) = 1
Solution (.solution) = [-0.999999, -0.999999, 4.0]
Number of iterations (.iter) = 28
Number of trust points (.p) = 4
Objective function value (.f) = 2.9377653725356374e-11
Number of function evaluations (.nf) = 28
Number of Jacobian evaluations (.nj) = 28
Index of outliers (.outliers) = [4]

Changing some options

RAFF has tunning options like precision of gradient stopping criteria and initial guess.

julia> output = raff(model, data, 3; initguess=[0.5,0.5,0.5], ε=1.0e-4)
** RAFFOutput **
Status (.status) = 1
Solution (.solution) = [-0.92378, -0.757813, -0.186976]
Number of iterations (.iter) = 28
Number of trust points (.p) = 5
Objective function value (.f) = 228.6473817033248
Number of function evaluations (.nf) = 28
Number of Jacobian evaluations (.nj) = 28
Index of outliers (.outliers) = Int64[]

RAFF is based on an optimization method. In this way, it is subject to stopping at stationary points that are not global minimizers. For this reason, heuristics were implemented to find global minimizers. Such heuristics depend on random number generation. So, if you want to run tests with more reliability this can be a useful strategy. To define in RAFF, say, 1000 different starting points, is enough to redefine the keyword argument MAXMS.

julia> output = raff(model, data, 3; MAXMS=1, initguess=[0.5,0.5,0.5], ε=1.0e-10)
** RAFFOutput **
Status (.status) = 1
Solution (.solution) = [-1.0, -1.0, 4.0]
Number of iterations (.iter) = 48
Number of trust points (.p) = 4
Objective function value (.f) = 4.536928392543292e-25
Number of function evaluations (.nf) = 48
Number of Jacobian evaluations (.nj) = 42
Index of outliers (.outliers) = [4]

In the above example, we have also changed the starting point for the method. Also, the stopping criterion was changed to $10^{-10}$, which means high accuracy when solving the subproblems. See RAFF API for all the possible options that can be used.

Parallel running

RAFF can be run in a parallel or distributed environment, using the Distributed package and function praff. Let's use praff to solve the same problem from the beginning. First, the Distributed package has to be loaded and the number of workers has to be added. It is also possible to add the address of other machines.

using Distributed

addprocs(3) # Add 3 worker processes

This step can be replaced if Julia is initialized with the -p option

julia -p 3

Now we have to load RAFF and the fit function in all workers:

@everywhere using RAFF

@everywhere function model(x, θ)
   θ[1] * x[1]^2 + θ[2]
end

@everywhere function gmodel!(g, x, θ)
    g[1] = x[1]^2
    g[2] = 1.0
end

then, we call praff to solve the problem (note that we do not need to send the A matrix to all workers, since it will be automatically sent by praff).

A=[-2.0  5.0;
  -1.5  3.25;
  -1.0  2.0 ;
  -0.5  1.25;
   0.0  1.0 ;
   0.5  2.55;
   1.0  2.0 ;
   1.5  3.25;
   2.0  5.0 ;];

n = 2

output = praff(model, gmodel!, A, n)
RAFFOutput(1, [1.0, 0.999996], 6, 8, 4.0205772365906425e-11, [6])

The true effectiveness of parallelism occurs when option MAXMS is set, which changes the number of random initial points that are tried for each subproblem solved. Better solutions can be achieved with higher values of MAXMS

n = 2

output = praff(model, gmodel!, A, n; MAXMS=1000)
RAFFOutput(1, [1.0, 1.0], 7, 8, 5.134133698545651e-13, [6])