Octave Wiki | RecentChanges

CodaTutorial

DaCodaAlFine - CodaTutorial

Octahedron, octopus, oculist, oct-file? -- What?

An oct-file is a dynamical extension of the Octave interpreter, in other words a shared library or shared object. The source files that make up an oct-file are written in C++ and therefore, conventionally, carry the extension cc.

Why would you ever want to write in C++ again, having the wonderful Octave environment at your service? Well, sometimes the performance of native Octave-code is not enough to attack a problem. Octave-code is interpreted, thus it is inherently slow when executed (especially if the problem cannot be vectorised). In such a case moving from m-code to compiled C++-code speeds up the execution by a factor of ten or more. The second group of problems, where reverting to low-level code is necessary, is when interfacing with foreign libraries (think of LAPACK) or routines written in languages other than C++, most notably Fortran-77.

Having convinced ourselves that we have to bite the bullet, we start with an easy CodaTutorial. This will teach any reader who is familiar with C++ how to write her first dynamically linked extension for Octave. Having guided the reader in her first steps, CodaAdvanced covers more advanced topics of writing extensions.

Tutorial

The tutorial unfurls slowly; we do not want to scare anybody by the complexity of the task ahead. Also look at the sources in OCT/src or OCT/src/DLD-FUNCTIONS, where OCT is the root directory of your Octave source-tree and consult the header-files and already existent dynamic extensions.

Problem Definition

Instead of giving an abstract treatise, we want to explain the whole business of dynamical extensions with the help of a simple, yet non-trivial problem. This example shall be analyzed in detail from beginning to end, and we shall elucidate some modern software construction principles on the way.

We shall implement a matrix power function. Given a square matrix A of integers, reals, or complex values and a non-negative integral exponent n, the function matpow(A, n) is to calculate An.

We have to watch out for boundary cases: an empty matrix A or a zero exponent n.

Note:
Octave already has a matrix power operator, the caret "^", which is more powerful than our example ever will be. This should not bother us.

High-Level Implementation

We postpone writing the C++-implementation until we are completely satisfied with our implementation in Octave. Having Octave, a rapid prototyping environment at hand, it would be stupid to throw away its possibilities.

Hobby Astronomer's Mirror Making Rule
It is faster to grind a 3 inch mirror and then a 5 inch mirror, than to grind a 5 inch mirror.

Now that the problem is exactly defined, we can start thinking about an implementation. Here is our naive first shot:

 function b = matpow(a, n)
     ## Return b = a^n for square matrix a, and non-negative, integral n.
     ## Note: matpow is an extremely naive implementation.
 
     b = eye(size(a));
     for i = 1:n
         b = b * a;
     endfor
 endfunction

Easy does the job! matpow looks like it does what we want, but how can we be sure? The solution lies in unit tests -- extremely useful (and necessary) when implementing code in C++.

Add the following snippet to the bottom of matpow:

 %!shared a
 %!test
 %!  a = [ 2.0, -3.0;
 %!       -1.0,  1.0];
 %!
 %!assert(matpow(a,0), diag([1,1]));
 %!assert(matpow(a,1), a);
 %!assert(matpow(a,2), a^2);
 %!assert(matpow(a,3), a^3);
 %!assert(matpow(a,4), a^4);
 %!assert(matpow(a,22), a^22);
 %!assert(matpow(a,23), a^23);

Running these tests on matpow gives us confidence!

 octave:1> test matpow
 PASSES 8 out of 8 tests

We now get more ambitious.

Our algorithm is really naive, and matrix multiplications are computationally expensive. Let us cut down on the number of multiplications. What is the minimum number of multiplications needed to compute the n-th power of A? Starting with A, we can square it, getting A². Again squaring is the fastest way to get to the fourth power. In general squaring our highest power lets us advance with fewest multiplications. This is the heart of our new algorithm.

Improved matrix power algorithm.
If the exponent n is a w-bit number, we can apply the binary expansion n = ∑(b[i] * 2i, i = 0..w-1), where the b[i] are either 0 or 1. In other words, we square A for every bit in the binary expansion of n, multiplying this value with the final result if bi = 1. Special care has to be taken for the cases where n = 0. See also [Golub:1996], section 11.2.5, page 569.

The new algorithm looks like this:

 function b = matpow(a, n)
     ## Return b = a^n for square matrix a, and non-negative, integral n.
 
     ## handle special case: n = 0
     if (n == 0)
         b = eye(size(a));
         return;
     endif
  
     ## general case: n >= 1
     p = a;                              # p holds the current square
     b = a;
     np = n - 1;
     while (np >= 1)
         if (is_even(np))                # is_even is defined below
             ## zero in the binary expansion
             np = np / 2;
         else
             ## one in the binary expansion
             np = (np - 1) / 2;
             b *= p;
         endif
         p *= p;
     endwhile
 endfunction
 
 function e = is_even(n)
     ## Return true if n is even, false otherwise.
 
     e = (rem(n, 2) == 0);
 endfunction

The new algorithm reduces the number of matrix multiplications from O(n), to O(log2(n)), which is a remarkable improvement for large matrices as matrix multiplication itself is an O(m³) process (m is the number of rows and columns of the matrix).

Running the test-suite again ensures that our code works as expected.

Our First Dynamic Extension

Ready for the jump to light speed? Wait, we have to feed the navigation computer first!

Feature Compatibility

In principle, the functions defined in oct-files must have the same capabilities as functions in m-files. In particular, the input and output arguments lists must be able to carry different numbers of arguments, which are moreover of different type, just like m-file functions can. This means that there must be a way of mapping a function from Octave code like

 function [array, real-scalar, integer] =
     func(complex-scalar, array, list, integer) 
     ## func frobs the snafu, returning all gromniz coefficients.
 
     -- actual function code goes here
 
 endfunction

to C++. To this end, Octave's core interface supplies the macro (in OCT/src/defun-int.h):

 octave_value_list
 DEFUN_DLD(function_name,
           const octave_value_list& input_argument_list,
           [ int number_of_output_arguments ],            -- optional argument
           const char* documentation_string)

We have decorated the macro arguments with their types. Note that the first argument is the name of the function to be defined.

Of course the macro has to be defined somewhere. The easiest way to pull in all necessary definitions is to include OCT/src/oct.h often installed as /usr/include/octave-VERSION/octave/oct.h. Now, our skeletal source file of the dynamic extension has the following shape:

 #include <octave/oct.h>
 
 DEFUN_DLD(func, args, nargout,
           "func frobs the snafu, returning the gromniz coefficients.")
 {
     /* -- actual function code goes here */
 
     return possibly_empty_list_of_values;
 }

Often functions do not depend on the actual number of output parameters. In this case the third argument to DEFUN_DLD can be omitted.

Tip - Naming Convention
DEFUN_DLD gives the user the freedom to choose any name for input_argument_list and number_of_output_arguments, but conventionally args and nargout are used (thus reminding of the parameter names in main, which are argc and argv).

Essential Types and Methods

Before we start to write a low-level implementation of matpow, we look at the most essential types and methods used to handle data that flows through the Octave-interpreter interface.

As has been said above, the arguments passed to dynamically loaded functions are bundled in an octave_value_list. Results returned from a function are also passed in an octave_value_list. The default constructor, octave_value_list, creates an empty list which, when used as return value, yields an Octave function returning "void". To access the elements of a list the following methods from class octave_value_list (defined in OCT/src/oct-obj.h) are useful:

 {
   public octave_value_list();
   public octave_value& operator()(int index);
   public const octave_value operator()(int index);
   public const int length();
 }

length returns the number of elements in the list. The overloaded parenthesis operator selects a single element from the list. The index base for index is 0 and not 1 as Fortran users might infer. The value returned by operator() is an octave_value, i.e., any type known to Octave, e.g. integers, real matrices, complex matrices, and so on.

Knowing how to access the arguments of unspecified type, the next thing we would like to do is get their values. Octave follows a uniform naming scheme: all functions that return an object of a certain type end in _type. Some of the more important of these methods (defined in OCT/src/ov-base.h) are:

 {
   public const int int_value();
   public const double double_value();
   public const Complex complex_value();
   public const Matrix matrix_value();
 }

Low-level implementation

Now, we are ready to implement matpow as a dynamically loaded extension to Octave:

 #include <octave/oct.h>
 
 static bool is_even(int n);
 
 DEFUN_DLD(matpow, args, ,
           "Return b = a^n for square matrix a, and non-negative, integral n.")
 {
     const int n = args(1).int_value();  
 
     if (n == 0)
         return octave_value(
                DiagMatrix(args(0).rows(), args(0).columns(), 1.0)
                );
 
     Matrix p(args(0).matrix_value());
     Matrix b(p);
     int np = n - 1;
     while (np >= 1)
     {
         if (is_even(np))
         {
             np = np / 2;
         }
         else
         {
             np = (np - 1) / 2;
             b = b * p;
         }
         p = p * p;
     }
 
     return octave_value(b);
 }
 
 bool
 is_even(int n)
 {
     return n % 2 == 0;
 }

Tip:
The Octave library overloads all elementary operations of scalars, (row-/column-) vectors and matrices. If in doubt as to whether a particular operation has been overloaded, simply try it. It takes less time than browsing (read: grepping through) the source code and the desired elementary operation is implemented in most cases.

We learn from the example that the C++ closely resembles the Octave function. This is due to the clever class structure of the Octave library.

Compiling

Now that we are done with the coding, we are ready to compile and link. The Octave distribution contains the script mkoctfile, which does exactly this for us. In the simplest case it is called with the C++-source as its only argument.

 $ ls
 Makefile  matpow.m  matpow.cc
 $ mkoctfile matpow.cc

Good! But while developing, we would like to compile several times, run our test suite and finally remove all files that can be re-generated from source. Enter: Makefile.

 # Makefile for Octave extension matpow
 
 .phony: all
 all: matpow.oct
 
 .phony: clean
 clean:
         rm -f matpow.oct *.o
 
 .phony: distclean
 distclean: clean
         rm -f *~ octave-core core
 
 .phony: test
 test: matpow.oct
         echo "test matpow.cc" | octave --silent
 
 %.oct: %.cc
         mkoctfile $<

Running

If the oct-file is in the LOADPATH, it will be loaded automatically -- either when requesting help on the function or when invoking it directly.

 $ ls
 Makefile  matpow.m  matpow.cc  matpow.o  matpow.oct  test_matpow.m
 $ octave -q
 octave:1> help matpow
 matpow is the dynamically-linked function from the file
 /home/cspiel/hsc/octave/src/matpow/matpow.oct
 
 Return b = a^n for square matrix a, and non-negative, integral n.
 
 octave:2> matpow([2, 0; 0, -1], 4)
 ans =
 
   16   0
    0   1
 
 octave:3> [Ctrl-D]

Spicing up matpow

The keen reader will have noticed that matpow is highly vulnerable to unexpected input: read on to see how this is handled! CodaAdvanced describes how to document our newly written function properly.

Input parameter checking

matpow will fail when called with an invalid set of parameters. In order to make it more robust, we examine the different parameters before the function is executed, making sure that they conform to our requirements.

Parameter tests should be done in the following order:

  1. Number of arguments
  2. Type of argument
  3. Size of argument (where applicable)
  4. Range of argument (where applicable and necessary)
  5. Inter-argument relations (if necessary)

Octave supplies the user with all necessary type- and size-tests as class member functions. The type-tests (defined in OCT/src/ov-base.h) share the common prefix is_. Here are a few:

 {
   public const bool is_real_scalar();
   public const bool is_complex_scalar();
   public const bool is_real_matrix();
   public const bool is_complex_matrix();
 }

Using these functions is seldomly necessary, however. Instead, we ask Octave to convert an input argument to the required type. If there is a problem (indicated by the error_state (OCT/src/error.h)), we abort, e.g.

  const int n = args(1).int_value();
  if (error_state) return octave_value_list();

To examine the sizes and shapes of different Octave objects, the following methods prove useful:

 {
   public const int rows();
   public const int columns();
   public const int length();
   public const bool is_empty();
   public const bool is_square(); 
 }

Remember that the methods available depend on the underlying type. For example, a ColumnVector only has a length (OCT/src/ov-list.h), whereas a Matrix has a number of rows and columns (OCT/src/ov-base-mat.h).

We have all the knowledge we need to modify matpow to include parameter checking:

    ...
    const int n = args(1).int_value();
    Matrix p(args(0).matrix_value());
    if (error_state) return retval;
    if (n < 0) {
        error("n must be positive");
        return retval;
    }
    if (!p.is_square()) {
        error("A must be square");
        return retval;
    }
    if (n == 0)
        return octave_value(
                            DiagMatrix(args(0).rows(), args(0).columns(), 1.0)
                            );
    Matrix b(p);
    int np = n - 1;   
    ...

Our final task is to update the tests and run them. Take the tests from matpow.m and append them to the C++ source inside a comment:

 /*
 %!shared a
 %!test
 %!  a = [ 2.0, -3.0;
 %!       -1.0,  1.0];
 ...
 %!assert(matpow(a,23), a^23);
 */

If you prefer to separate your tests from your code, you can put them in a file called 'matpow' instead (remember to update the Makefile).

Now, also add the following tests:

 %!fail("matpow([1,1; 1,1])", "usage:");
 %!fail("matpow()", "usage:");
 %!fail("matpow([1,1; 1 1], 2, 1)", "usage:");
 %!fail("matpow([1 2 3; 4 5 6], 1)");
 %!fail("matpow(a, -1)");


 $ octave -q
 octave:1> test matpow.cc
 PASSES 13 out of 13 tests

Texinfo documentation strings

Our much improved matpow still carries around the poor documentation string:

 "Return b = a^n for square matrix a, and non-negative, integral n."

Let us improve on this! Octave supports documentation strings -- docstrings for short -- in Texinfo format. The effect on the online documentation will be small, but the appearance of printed documentation will be greatly improved.

The fundamental building block of Texinfo documentation strings is the Texinfo-macro @deftypefn, which takes two arguments: The class the function is in, and the function's signature. For functions defined with DEFUN_DLD, the class is Loadable Function. A skeletal Texinfo docstring therefore looks like this:

 -*- texinfo -*- 
 @deftypefn{Loadable Function} {return_values =} function_name(arguments)
 
 @cindex index term 
 
 Help text in Texinfo format...
 @end deftypefn

Texinfo has several macros which control the markup. In the group of format changing commands, we note @var{variable_name}, and @code{code_piece}. The Texinfo format has been designed to generate output for online viewing with text-terminals as well as generating high-quality printed output. To these ends, Texinfo has commands which control the diversion of parts of the document into a particular output processor. Two formats are of importance "info" and "TeX". The first is selected with

 @ifinfo
 text for info only
 @end ifinfo

the latter with

 @iftex
 @tex
 text for TeX only
 @end tex
 @end iftex

If no specific output processor is chosen, by default the text goes into both (or, more precisely: all) output processors. Usually, neither @ifinfo, nor @iftex appear alone, but always in pairs, as the same information must be conveyed in every output format.

Here is a sample docstring for matpow.cc:

 -*- texinfo -*-
 @deftypefn{Loadable Function} {@var{b} =} matpow(@var{a}, @var{n})
 
 @cindex matrix power
 
 Return matrix @var{a} raised to the @var{n}-th power.  Matrix @var{a} is
 a square matrix and @var{n} a non-negative integer.
 @iftex 
 @tex 
 $n = 0$
 @end tex
 @end iftex
 @ifinfo 
 n = 0
 @end ifinfo
 is explicitly allowed, returning a unit-matrix of the same size as
 @var{a}.
 
 Note: @code{matpow}  duplicates part of the functionality of the built-in
 exponentiation operator
 @iftex
 @tex
 ``$\wedge$''.
 @end tex
 @end iftex
 @ifinfo
 @code{^}.
 @end ifinfo
 
 
 Example:
 @example 
 matpow([1, 2; 3, 4], 4)
 @end example
 @noindent
 returns
 @example
 ans =
 
   199  290
   435  634
 @end example
 
 The algorithm takes
 @iftex
 @tex
 $\lfloor \log_{2}(n) \rfloor$
 @end tex
 @end iftex
 @ifinfo
 floor(log2(n))
 @end ifinfo
 matrix multiplications to compute the result.
 @end deftypefn

For further information about Texinfo consult the Texinfo documentation. For TeX-beginners we recommend "The Not So Short Introduction to LaTeX" by Tobias Oetiker et. al.

One thing we held back is the true appearance of a Texinfo docstring -- mainly because it looks so ugly. The C++-language imposes the constraint that the docstring must be a string-constant. Moreover, because DEFUN_DLD is a macro, every line-end has to be escaped with a backslash. The backslash does not insert any whitespace and TeX separates paragraphs with empty lines, so that we have to put in new-lines as line-separators. Thus, the Texinfo docstring in source form has each line end decorated with "\n\".

 DEFUN_DLD(matpow, args, ,
           "-*- texinfo -*-\n\
 @deftypefn{Loadable Function} {@var{b}} = matpow(@var{a}, @var{n})\n\
 \n\
 @cindex matrix power\n\
 \n\
 Return matrix @var{a} raised to the @var{n}-th power.  Matrix @var{a} is\n\
 a square matrix, and @var{n} a non-negative integer.\n\
 ...")

At least the formatted versions look much better. The info-version (which will be used in Octave's online help) looks like this:

 - Loadable Function: B = matpow(A, N)
     Return matrix A raised to the N-th power.  Matrix A is a square
     matrix and N a non-negative integer.  n = 0 is explicitly allowed,
     returning a unit-matrix of the same size as A.
 
     Note: `matpow'  duplicates part of the functionality of the
     built-in exponentiation operator `^'.
 
     Example:
          matpow([1, 2; 3, 4], 4)
 
     returns
          ans =
            199  290
            435  634
 
     The algorithm takes floor(log2(n)) matrix multiplications to
     compute the result.

DaCodaAlFine - CodaTutorial

Octave Wiki | RecentChanges
Edit text of this page | View other revisions
Last edited 2005-11-15 3:56 UTC (diff)
Search:

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation