...making Linux just a little more fun!
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 endfunctionListing 1: lg_factorial1.m
So, a function definition in Octave is:
function name body endfunctionFunctions 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/:"
~/.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
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 endfunctionListing 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 endfunctionListing 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 endfunctionListing 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 endfunctionListing 5: lg_factorial5.m
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 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>
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.
## 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 endfunctionListing 8: lg_factorial.m
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 to and if I get
enough I will write them up in an article.
Barry O'Donovan graduated from the National University of Ireland, Galway
with a B.Sc. (Hons) in computer science and mathematics. He is currently
completing a Ph.D. in computer science with the Information Hiding Laboratory, University
College Dublin, Ireland in the area of audio watermarking.
Barry has been using Linux since 1997 and his current flavor of choice
is Fedora Core. He is a member of the Irish
Linux Users Group. Whenever he's not doing his Ph.D. he can usually be
found supporting his finances by doing some work for Open Hosting, in the pub with friends or running in the local
park.