BarryODonovan.com :: GNU Octave - Functions and Scripts

GNU Octave - Functions and Scripts

Originally published in the LinuxGazette.net, March 2005, Issue 112.

Introduction

In this, the second article on GNU Octave, I will build on the basics that were covered in issue #109 by introducing functions and scripts through a number of examples. Obtaining and installing GNU Octave as well as sources for official documentation is discussed in the first article and you should refer to it for more information.

Functions in Octave

Just like any other programming language, Octave has full support for creating functions. Functions are an essential tool that allow large problems to be broken down into a number of smaller tasks. A function should perform a specific task and it should perform it well. These criteria are very important. The more specific the task that a function performs the more reusable it will be; although you may be writing it to help solve your current problem, if it is well defined then it could also be used in many future problems. By "perform it well", I mean that the function should give the correct answer for valid input while reporting an error for invalid input; it should be a veritable "black-box" - once written and tested it should be of sufficient quality that you can trust it without having to recheck it in all future problems.

One extremely simply mathematical function that Octave lacks and which I required recently was the factorial function. This function is defined as:

n! = n(n-1)...2·1

So, for example,

5! = 5·4·3·2·1 =
120

. There are a number of algorithms for implementing this function and I will cover some of them to give a good introduction to functions in Octave. Let's begin with an iteration solution:

function answer = lg_factorial1( n )

    answer = 1;
    
    for i = 2:n
        answer = answer * i;
    endfor

endfunction

Listing 1: lg_factorial1.m

So, a function definition in Octave is:

function name
    body
endfunction

Functions should be saved on their own in a text file named the same as the function itself with the extension .m. So the above file would be saved as lg_factorial1.m. When you now try and call the function lg_factorial1(), Octave will search the list of directories specified by the built-in variable LOADPATH for files ending in ‘.m’ that have the same base name as the function name. If you want to create a repository of functions on your computer and would like to have LOADPATH include that directory automatically, you can add the line:
LOADPATH = "/path/to/your/files/:"
to the Octave configuration file ~/.octaverc. Octave also checks the path specified in the built-in variable DEFAULT_LOADPATH which includes the current working directory by default.

A function can take any number of arguments as a comma separated list in parentheses after the function name. Multiple return values can also be defined:

function [retval1, retval2, etc] = name( arg1, arg2, etc )
    body
endfunction

There are two additional rules in the mathematical definition of the factorial function:

0! = 1
n >= 0

Let's incorporate these rules into our function definition:

function answer = lg_factorial2( n )

    if( n < 0 )
        error( "there is no definition for negative factorials" );
    endif
    
    answer = 1;
    
    if( n == 0 )
        return;
    else
        for i = 2:n
            answer = answer * i;
        endfor
    endif

endfunction

Listing 2: lg_factorial2.m The function first tests to ensure that the input is valid (non-negative). If it is not it throws an error using the error() built-in function. As well as printing the error message, it also prints a traceback of of all the functions leading to the error. This is very useful for programmers who are solving complex problems with many functions as they can narrow done the offending problem very quickly.

If the input is valid then we test for the zero case and we use the return command to end the function if it is true. Unlike many other programming languages, Octave's return statement does not take any arguments, so it is essential that the return value(s) are set before a return call is encountered as I have done with the answer = 1 statementabove.

Now, what happens if lg_factorial2() is called without any arguments?

octave:1> lg_factorial2()
error: `n' undefined near line 3 column 9
error: evaluating binary operator `<' near line 3, column 11
error: if: error evaluating conditional expression
error: evaluating if command near line 3, column 5
error: called from `lg_factorial2' in file `/home/user/lg/lg_factorial2.m'
octave:1>

We can also check to ensure that a valid number of arguments have been passed to the function using the built-in variable nargin. This variable is automatically initialised to the number of arguments passed. While we're at it, we should also try and ensure that a valid data type is passed. Unfortunately, Octave has no isinteger() function but we can check that it is a real number and not a vector. If a real non-integer is passed it will be automatically rounded down by the range operator (2:n).

function answer = lg_factorial3( n )

    if( nargin != 1 )
        usage( "factorial( n )" );
    elseif( !isscalar( n ) ||  !isreal( n ) )
        error( "n must be a positive integer value" );
    elseif( n < 0 )
        error( "there is no definition for negative factorials" );
    endif
    
    if( n == 0 )
        answer = 1;
        return;
    else
        answer = prod( 1:n );
    endif

endfunction

Listing 3: lg_factorial3.m

This example introduces the usage() built-in function. It is very similar to error() in that it prints a traceback of the functions leading up to the call to help debugging but instead of printing "error: ...", it prints "usage: ...". The isscalar() and isreal() functions check that the given argument is a scalar value (as opposed to a vector, a string, etc) and a real number (as opposed to a complex number) respectively. The ! before them inverts the test so that it reads as if n is not a scalar or if n is not a real number then raise an error.

You also have noticed that I have changed the code for calculating the factorial. I now use the built-in function prod() which calculates the product of the elements in a given vector. In this case, the given vector is the range 1:n.

In the previous article, I mentioned that Octave comes with full documentation which is accessible through the Linux command info. There is also a built-in help function for every command. For example:

octave:2> help prod
prod is a built-in function

 -- Built-in Function:  prod (X, DIM)
     Product of elements along dimension DIM.  If DIM is omitted, it
     defaults to 1 (column-wise products).
...
octave:3>

A "well written" function should allow the use of the help function in a similar manner. The help text is taken as the first block of non-copyright (see below) comments from a function file:

## usage: answer = lg_factorial4( n )
##
## Returns the factorial of n (n!). n should be a positive
## integer or 0.

function answer = lg_factorial4( n )

    if( nargin != 1 )
        usage( "factorial( n )" );
    elseif( !isscalar( n ) ||  !isreal( n ) )
        error( "n must be a positive integer value" );
    elseif( n < 0 )
        error( "there is no definition for negative factorials" );
    endif
    
    if( n == 0 )
        answer = 1;
        return;
    else
        answer = prod( 1:n );
    endif

endfunction

Listing 4: lg_factorial4.m

And calling the help function will now yield:

octave:3> help lg_factorial4
lg_factorial4 is the user-defined function from the file
/home/user/lg/lg_factorial4.m

usage: answer = lg_factorial4( n )

Returns the factorial of n (n!). n should be a positive
integer or 0.

...
octave:4>

Another common algorithm for calculating the factorial of a number is to use a recursive function (a recursive function is one that calls itself). This can be implemented in Octave (ignoring the error checking for now) as:

function answer = lg_factorial5( n )

    if( n == 0 )
        answer = 1;
        return;
    else
        answer = n * lg_factorial5( n -1 );
    endif

endfunction

Listing 5: lg_factorial5.m

Testing Function Efficiency

Academic use of Octave often tends to be in the form of large and time-consuming simulations. As such it is important to know how to both write quick and efficient code as well as testing the efficiency of code that you do write.

Let's begin by comparing the three algorithms above for calculating the factorial of a number. Octave has two functions for starting and stopping a "wall-clock timer":

octave:4> tic(); sleep( 10 ); toc()
ans = 10.0
octave:5>

To make the comparison fair, the only checking we will perform is if n == 0, otherwise we will assume it is a positive integer. We will use the program in listing 5 above, the iterative version shown in lg_factorial6.m and the prod() version shown in lg_factorial7.m. The commands executed for the comparison and the times found can be seen here:

octave:5> tic(); for i=1:100 for n=1:100 lg_factorial5( n ); end; end; toc()
ans = 23.5154919996858
octave:6> tic(); for i=1:100 for n=1:100 lg_factorial6( n ); end; end; toc()
ans = 3.06905199959874
octave:7> tic(); for i=1:100 for n=1:100 lg_factorial7( n ); end; end; toc()
ans = 0.537685000337660
octave:8>

Firstly, the recursive function took the longest by far and this was to be expected. Just like with most other programming languages, calling functions (even recursively) requires a lot of overhead. The result that may surprise you is that the use of the prod() function is about six times faster than iteration. This is because, just like Matlab, Octave is quite slow at looping and it should be avoided wherever possible.

As I have already stated, well written and documented functions are important, especially if you wish to share your code. Some may think that the effect of checking for valid input may significantly slow a function's execution. Let's compare the iteration algorithm in lg_factorial6.m against the same algorithm but with full error checking in lg_factorial8.m:

octave:8> tic(); for i=1:100000 lg_factorial6( 10 ); end; toc()
ans = 9.44780500046909
octave:9> tic(); for i=1:100000 lg_factorial8( 10 ); end; toc()
ans = 17.9307480007410
octave:10>

There is clearly a difference but when you consider that each function is called 100,000 times, the added time for the error checking is only 0.000085 seconds per function call.

Another important tip to remember when writing Octave functions is to avoid resizing matrices unnecessarily. The manual itself states that if you are “building a single result matrix from a series of calculations, set the size of the result matrix first, then insert values into it.” The following is a graphic illustration at just how slow this can be:

octave:10> clear a; tic(); for i=1:100000 a(i) = i; end; toc()
ans = 52.4
octave:11> a = [1]; tic(); for i=2:100000 a = [a i]; end; toc()
ans = 362.247
octave:12> a=zeros(100000,1); tic(); for i=1:100000 a(i) = i; end; toc()
ans = 1.42
octave:13>

Octave Scripts

Octave script files are simply a sequence of Octave commands which are evaluated as if you had typed them at the Octave prompt yourself. They are useful for setting up simulations where you wish to vary certain parameters for each run without having to re-type the commands each time, for sequences of commands that do not logically belong in a function and for automating certain tasks.

Octave scripts are placed in a file with the .m extension in the same way that functions are but a script file must not begin with the function keyword. Note that variables defined in a script share the same namespace (or scope) as variables defined at the Octave prompt.

The following is a simple example of an Octave script which calculates the time taken to create an array of integers of size specified by the user:

# An example Octave script 

len = input( "What size array do you wish to use for the evaluation: " );

clear a; 
tic(); 
for i=1:len
    a(i) = i; 
endfor 
time1 = toc();

a = [1]; 
tic(); 
for i=2:len 
    a = [a i]; 
endfor
time2 = toc();

a=zeros( len, 1 ); 
tic(); 
for i=1:len 
    a(i) = i; 
endfor
time3 = toc();

printf( "The time taken for method 1 was %.4f seconds\n", time1 );
printf( "The time taken for method 2 was %.4f seconds\n", time2 );
printf( "The time taken for method 3 was %.4f seconds\n", time3 );

Listing 6: lg_calc_time.m

And once it is saved in an appropriately named file, it can be executed as follows:

octave:13> lg_calc_time
What size array do you wish to use for the evaluation: 20000
The time taken for method 1 was 1.3509 seconds
The time taken for method 2 was 12.8984 seconds
The time taken for method 3 was 0.2850 seconds
octave:14>

Executable Octave Scripts

It is also possible to have executable Octave script files like we have executable Bash scripts. This is an excellent feature of Octave that can be very useful for large problems where a mixture of programs, tools or applications may be needed to solve the problem.

The example in listing 6 above can be easily converted to an executable Octave program by adding a single line to the beginning:

#! /usr/bin/octave -qf

# An example executable Octave script 

len = input( "What size array do you wish to use for the evaluation: " );
... ... ...              ... ... ...              ... ... ...
printf( "The time taken for method 3 was %.4f seconds\n", time3 );

Listing 7: lg_calc_time.sh

Just make sure that the path to the Octave program is correct for your system (/usr/bin/octave), make the script executable and run it as follows:

[barry@hiscomputer lg]$ chmod u+x lg_calc_time.sh
[barry@hiscomputer lg]$ ./lg_calc_time.sh
What size array do you wish to use for the evaluation: 20000
The time taken for method 1 was 1.3959 seconds
The time taken for method 2 was 13.0201 seconds
The time taken for method 3 was 0.2800 seconds
[barry@hiscomputer lg]$

If you call the script with command line arguments, then they will be available through the built-in variable argv and the number of arguments will be contained in the variable nargin which we have already seen. A simple example of this can be seen in lg_calc_time2.sh which is the same as the last example except it reads the size of the array from the command line.

Conventional Headers and Copyright Notices

As I mentioned in the previous article, Octave is not as complete as Matlab for specialised functions and the Octave developers welcome new additions which are well written and robust. There are conventions for writing function files which include author information, copyright restrictions, date of creation, version, etc. Functions for inclusion with Octave will have to be distributed with an appropriate open source license (further information can be found in appendix A of the official documentation). With proper headers and a copyright notices, my factorial function would look like:

## Copyright (C) 2005 Barry O'Donovan
##
## 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, or (at your option) any later version.
##
## Octave 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.
##
## You should have received a copy of the GNU General Public
## License along with Octave; see the file COPYING.  If not,
## write to the Free Software Foundation, 59 Temple Place -
## Suite 330, Boston, MA 02111-1307, USA.

## usage: answer = lg_factorial( n )
##
## Returns the factorial of n (n!). n should be a positive
## integer or 0.

## Author: Barry O'Donovan <barry@ihl.ucd.ie>
## Maintainer: Barry O'Donovan <barry@ihl.ucd.ie>
## Created: February 2005
## Version: 0.1
## Keywords: factorial

function answer = lg_factorial( n )

    if( nargin != 1 )
        usage( "factorial( n )" );
    elseif( !isscalar( n ) ||  !isreal( n ) )
        error( "n must be a positive integer value" );
    elseif( n < 0 )
        error( "there is no definition for negative factorials" );
    endif
    
    if( n == 0 )
        answer = 1;
        return;
    else
        answer = prod( 1:n );
    endif

endfunction

Listing 8: lg_factorial.m

Final Words

There are many tricks and techniques for writing efficient Octave code. Examples include imaginative uses of the reshape() function for operations between elements of the same matrix so as to avoid iteration. I would welcome all such tips by e-mail and if I get enough I will write them up in an article.

BarryODonovan.com :: Publications