Document optimization.

Jürgen Knödlseder
1 parent 13f15a9c
Showing 7 changed files with 42 additions and 11 deletions
doc/source/user_manual/modules/linalg.rst
doc/source/user_manual/modules/opt.rst
doc/source/user_manual/modules/uml_opt.png
examples/cpp/optimize/optimize.cpp
include/GOptimizerPars.hpp
pyext/GOptimizerPars.i
src/opt/GOptimizerPars.cpp
-.. _sec_linalg:

Linear algebra
--------------

Overview
~~~~~~~~

The linear algebra module provides classes for vector and matrix 
manipulation. The ``GVector`` class implements a double precision
floating point vector, the classes ``GMatrix``, ``GMatrixSymmetric``
and ``GMatrixSparse`` implement a matrix of double precision
floating point values for a generic matrix, a symmetric matrix and
a sparse matrix, respectively. All matrix classes derive from the
abstract ``GMatrixBase`` base class that provides common functionalities
to all matrix classes.

:ref:`fig_uml_linalg` present an overview over the C++ classes of
the linear algebra module and their relations.

.. _fig_uml_linalg:

.. figure:: uml_linalg.png
   :width: 100%

   Linear algebra module

The classes ``GSparseNumeric`` and ``GSparseSymbolic`` are classes
that are used by ``GMatrixSparse`` but these classes are not
exposed to the outside world (i.e. the class definitions are not
part of the GammaLib interface).


Matrix Storage classes
~~~~~~~~~~~~~~~~~~~~~~

Three matrix storage classes are implemented:

* ``GMatrix`` implements a plain matrix storage class which stores
  all elements of the matrix in memory;
* ``GMatrixSymmetric`` implements a symmetric matrix storage class
  which stores only a triangle of matrix elements, imposing thus
  strict matrix symmetry;
* ``GMatrixSparse`` implements a sparse matrix storage class which
  stores only non-zero matrix elements in a column-wise organisation.

Storage class conversion constructors exist for all three classes
to transform one storage class into another::

    1  GMatrix          plain(10,10);
    2  GMatrixSymmetric symmetric(plain);
    3  GMatrixSparse    sparse(symmetric);

In the above example, a plain matrix is instantiated in line 1, the
plain matrix is converted into a symmetric matrix in line 2, and the
symmetric matrix is converted into a sparse matrix in line 3.
Additional complementary storage class conversion constructors exist,
but conversion to a symmetric matrix is of course only possible if the
matrix is indeed symmetric.

Matrix elements are accessed using the ``operator()``. Filling of
sparse matrix elements using this operator is possible, although this
can be time consuming due to internal memory management. In general,
each fill of a new non-zero element needs to shift all elements
that are located after that element in memory.

To reduce the memory management overhead in the filling of a sparse 
matrix, methods have been implemented that allow to fill a
matrix column wise::

    1  GMatrixSparse sparse(10,5);
    2  GVector       column(10);
    3  column[0] = 1.0;
    4  column[1] = 2.0;
    5  column[5] = 8.0;
    6  sparse.column(0, column);
    7  sparse.add_to_column(0, column);

Line 1 allocates a sparse matrix with 10 rows and 5 columns, line 2
instantiates a vector with 10 elements. In lines 3-5, 3 elements of
the vector are set to specific values, all other elements will default
to 0. In line 6, the elements of the vector are set as the elements
of the first matrix column (column 0). Line 7 differs from line 6 in
that the elements are now not set but added to the existing matrix
elements.

To further reduce the memory management overhead for the column-wise
fill of a sparse matrix, a "fill-stack" has been implemented. The
"fill-stack" is a buffer that implements a queue for columns that are
to be set or added to the matrix. The columns will be stored in this
"fill-stack" in the order they are provided, and only once the 
"fill-stack" is full, or upon request, the "fill-stack" will be flushed
into memory. The "fill-stack" is used as follows::

    sparse.stack_init(size, entries);
    ...
    sparse.column(0, column);
    ...
    sparse.stack_flush();
    ...
    sparse.stack_destroy();

The ``stack_init(size, entries)`` method initialises the "fill-stack",
where ``size`` is the size of the allocated memory buffer and ``entries``
is the maximum number of columns that will be held by the buffer.
If ``size`` is set to 0 (the default value), a default ``size`` value of
512 is used. If ``entries`` is set to 0 (the default value), the number of
matrix columns is taken as default ``entries`` value. Note that a too large
number of elements will produce some overhead due to "fill-stack"
management, hence ``entries`` should not exceeed a value of the order of
10-100.

The ``stack_flush()`` method flushes the stack, which is mandatory
before any usage of the matrix. Note that the "fill-stack" **is not
inserted automatically** before any matrix operation, hence manual stack
flushing is needed to make all filled matrix elements available for usage.
The ``stack_destroy()`` method will flush the stack and free all stack
elements. This method should be called once no filling is required anymore.
If ``stack_destroy()`` is called immediately after filling, no call to 
``stack_flush()`` is needed as the ``stack_destroy()`` method flushes the
stack before destroying it. The matrix stack is also destroyed by the
sparse matrix destructor, hence manual stack destruction is not
mandatory.
 \ No newline at end of file
+.. _sec_linalg:

Linear algebra
--------------

Overview
~~~~~~~~

The linear algebra module provides classes for vector and matrix 
manipulation. The ``GVector`` class implements a double precision
floating point vector, the classes ``GMatrix``, ``GMatrixSymmetric``
and ``GMatrixSparse`` implement a matrix of double precision
floating point values for a generic matrix, a symmetric matrix and
a sparse matrix, respectively. All matrix classes derive from the
abstract ``GMatrixBase`` base class that provides common functionalities
to all matrix classes.

:ref:`fig_uml_linalg` presents an overview over the C++ classes of
the linear algebra module and their relations.

.. _fig_uml_linalg:

.. figure:: uml_linalg.png
   :width: 100%

   Linear algebra module

The classes ``GSparseNumeric`` and ``GSparseSymbolic`` are classes
that are used by ``GMatrixSparse`` but these classes are not
exposed to the outside world (i.e. the class definitions are not
part of the GammaLib interface).


Matrix Storage classes
~~~~~~~~~~~~~~~~~~~~~~

Three matrix storage classes are implemented:

* ``GMatrix`` implements a plain matrix storage class which stores
  all elements of the matrix in memory;
* ``GMatrixSymmetric`` implements a symmetric matrix storage class
  which stores only a triangle of matrix elements, imposing thus
  strict matrix symmetry;
* ``GMatrixSparse`` implements a sparse matrix storage class which
  stores only non-zero matrix elements in a column-wise organisation.

Storage class conversion constructors exist for all three classes
to transform one storage class into another::

    1  GMatrix          plain(10,10);
    2  GMatrixSymmetric symmetric(plain);
    3  GMatrixSparse    sparse(symmetric);

In the above example, a plain matrix is instantiated in line 1, the
plain matrix is converted into a symmetric matrix in line 2, and the
symmetric matrix is converted into a sparse matrix in line 3.
Additional complementary storage class conversion constructors exist,
but conversion to a symmetric matrix is of course only possible if the
matrix is indeed symmetric.

Matrix elements are accessed using the ``operator()``. Filling of
sparse matrix elements using this operator is possible, although this
can be time consuming due to internal memory management. In general,
each fill of a new non-zero element needs to shift all elements
that are located after that element in memory.

To reduce the memory management overhead in the filling of a sparse 
matrix, methods have been implemented that allow to fill a
matrix column wise::

    1  GMatrixSparse sparse(10,5);
    2  GVector       column(10);
    3  column[0] = 1.0;
    4  column[1] = 2.0;
    5  column[5] = 8.0;
    6  sparse.column(0, column);
    7  sparse.add_to_column(0, column);

Line 1 allocates a sparse matrix with 10 rows and 5 columns, line 2
instantiates a vector with 10 elements. In lines 3-5, 3 elements of
the vector are set to specific values, all other elements will default
to 0. In line 6, the elements of the vector are set as the elements
of the first matrix column (column 0). Line 7 differs from line 6 in
that the elements are now not set but added to the existing matrix
elements.

To further reduce the memory management overhead for the column-wise
fill of a sparse matrix, a "fill-stack" has been implemented. The
"fill-stack" is a buffer that implements a queue for columns that are
to be set or added to the matrix. The columns will be stored in this
"fill-stack" in the order they are provided, and only once the 
"fill-stack" is full, or upon request, the "fill-stack" will be flushed
into memory. The "fill-stack" is used as follows::

    sparse.stack_init(size, entries);
    ...
    sparse.column(0, column);
    ...
    sparse.stack_flush();
    ...
    sparse.stack_destroy();

The ``stack_init(size, entries)`` method initialises the "fill-stack",
where ``size`` is the size of the allocated memory buffer and ``entries``
is the maximum number of columns that will be held by the buffer.
If ``size`` is set to 0 (the default value), a default ``size`` value of
512 is used. If ``entries`` is set to 0 (the default value), the number of
matrix columns is taken as default ``entries`` value. Note that a too large
number of elements will produce some overhead due to "fill-stack"
management, hence ``entries`` should not exceeed a value of the order of
10-100.

The ``stack_flush()`` method flushes the stack, which is mandatory
before any usage of the matrix. Note that the "fill-stack" **is not
inserted automatically** before any matrix operation, hence manual stack
flushing is needed to make all filled matrix elements available for usage.
The ``stack_destroy()`` method will flush the stack and free all stack
elements. This method should be called once no filling is required anymore.
If ``stack_destroy()`` is called immediately after filling, no call to 
``stack_flush()`` is needed as the ``stack_destroy()`` method flushes the
stack before destroying it. The matrix stack is also destroyed by the
sparse matrix destructor, hence manual stack destruction is not
mandatory.
 \ No newline at end of file
-.. _sec_opt:

Optimizers
----------
 \ No newline at end of file
+.. _sec_opt:

Optimizers
----------

Overview
~~~~~~~~

The optimizer module provides classes for function optimization.
The abstract ``GOptimizerFunction`` base class defines the interface
for the function that should be optimized. The ``GOptimizerPar``
class defines a parameter of the function, and the ``GOptimizerPars``
container class collects all function parameters. The optimizer
is represented by the abstract ``GOptimizer`` base class. So far,
the only optimizer algorithm implemented in GammaLib is the
Levenberg-Marquardt algorithm, implemented by the class
``GOptimizerLM``.


:ref:`fig_uml_opt` presents an overview over the C++ classes of
the optimizer module and their relations.

.. _fig_uml_opt:

.. figure:: uml_opt.png
   :width: 100%

   Optimizer module

The model parameter class ``GModelPar``, as a specific implementation
of the optimizer parameter class, derives from ``GOptimizerPar``.
One implementation of an optimizer function is the
``GObservations::likelihood`` class that is used for maximum
likelihood fitting within GammaLib.


Optimizing a function
~~~~~~~~~~~~~~~~~~~~~

This example illustrates how the minimum of a quadratic function of 
the form :math:`f(x)=2ax^2+bx+c` can be determined using the optimizer.
Please see ``examples/cpp/optimize`` for the source code.

First, the function is implemented as a class derived from
``GOptimizerFunction``::

    1  class function : public GOptimizerFunction {
    2  public:
    3      function(void) : m_value(0), m_gradient(1), m_covar(1,1) {}
    4      void           eval(const GOptimizerPars& pars);
    5      double         value(void) { return m_value; }
    6      GVector*       gradient(void) { return &m_gradient; }
    7      GMatrixSparse* covar(void) { return &m_covar; }
    8  protected:
    9      double        m_value;    //!< Function value
   10      GVector       m_gradient; //!< Function gradient vector
   11      GMatrixSparse m_covar;    //!< Covariance matrix
   12  };
   13  void function::eval(const GOptimizerPars& pars)
   14  {
   15      const double a =  2.0;
   16      const double b = -4.0;
   17      const double c =  2.0;
   18      double x       = pars[0]->value();
   19      m_value        = a*x*x + b*x + c;
   20      m_gradient[0]  = 2.0*a*x + b;
   21      m_covar(0,0)   = m_gradient[0] * m_gradient[0];
   22  };

Lines 1-12 define the ``function`` class that requires implementing the
``eval``, ``value``, ``gradient`` and ``covar`` methods. The ``eval``
method, implement in lines 13-22, performs the computation of the
function value, the gradient and the products of the gradients.

The optimization is then done using the following code::

    1  GOptimizerLM opt;
    2  function fct;
    3  GOptimizerPars pars(1);
    4  pars[0]->value(1.5);
    5  opt.optimize(fct, pars);
    6  std::cout << "Function value .....: " << fct.value() << std::endl;
    7  std::cout << "Parameter value ....: " << pars[0]->value() << std::endl;

Line 1 allocates an Levenberg-Marquardt optimizer, line 2 creates an 
instance of the function to optimize. In line 3, a parameter container
with a single parameter is allocated, and the value of the single parameter
is set to 1.5 in line 4. In line 5, the optimizer is called, and the 
resulting function value and best fitted parameter is logged in to console 
in lines 6-7.
 \ No newline at end of file
@@ -47,12 +47,13 @@ protected:
 };
 void function::eval(const GOptimizerPars& pars)
 {
-    const double a = 2.0;
-    const double b = 5.0;
-    const double c = 1.0;
-    double x       = pars.par(0).value();
+    const double a =  2.0;
+    const double b = -4.0;
+    const double c =  2.0;
+    double x       = pars[0]->value();
     m_value        = a*x*x + b*x + c;
-    m_gradient[0]  = 2*a*x + b;
+    m_gradient[0]  = 2.0*a*x + b;
+    m_covar(0,0)   = m_gradient[0] * m_gradient[0];
 };
  
  
@@ -63,22 +64,27 @@ void function::eval(const GOptimizerPars&amp; pars)
  ***************************************************************************/
 int main(void) {
  
+    // Allocate logger
+    GLog log;
+    log.cout(true);
+    
     // Allocate optimizer
+    //GOptimizerLM opt(log);
     GOptimizerLM opt;
  
     // Allocate function
-    function fct();
+    function fct;
  
     // Allocate parameters and set initial value
     GOptimizerPars pars(1);
-    pars[0].value(3.0);
+    pars[0]->value(1.5);
  
     // Optimize parameters
     opt.optimize(fct, pars);
  
     // Print derivatives
-    std::cout << "Function value: " << fct.value() << std::endl;
-    std::cout << "Parameter ....: " << pars[0].value() << std::endl;
+    std::cout << "Function value .....: " << fct.value() << std::endl;
+    std::cout << "Parameter value ....: " << pars[0]->value() << std::endl;
  
     // Exit
     return 0;
@@ -59,6 +59,7 @@ class GOptimizerPars : public GContainer {
 public:
     // Constructors and destructors
     GOptimizerPars(void);
+    explicit GOptimizerPars(const int& number);
     GOptimizerPars(const GOptimizerPars& pars);
     virtual ~GOptimizerPars(void);
  
@@ -38,6 +38,7 @@ class GOptimizerPars : public GContainer {
 public:
     // Constructors and destructors
     GOptimizerPars(void);
+    explicit GOptimizerPars(const int& number);
     GOptimizerPars(const GOptimizerPars& pars);
     virtual ~GOptimizerPars(void);
  
@@ -64,6 +64,29 @@ GOptimizerPars::GOptimizerPars(void)
  
  
 /***********************************************************************//**
+ * @brief Parameter constructor
+ *
+ * @param[in] number Number of parameters to allocate.
+ *
+ * Allocates @p number parameters in the parameter container.
+ ***************************************************************************/
+GOptimizerPars::GOptimizerPars(const int& number)
+{
+    // Initialise members
+    init_members();
+
+    // Allocate parameters
+    for (int i = 0; i < number; ++i) {
+        m_alloc.push_back(true);
+        m_pars.push_back(new GOptimizerPar);
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Copy constructor
  *
  * @param[in] pars Optimizer parameters.