LINUX GAZETTE
[ Prev ][ Table of Contents ][ Front Page ][ Talkback ][ FAQ ][ Next ]

"Linux Gazette...making Linux just a little more fun!"

Numerical Workbenches, part II

By Christoph Spiel

In Part I , we looked at the most basic operations of the numerical workbenches GNU/Octave 2.1.34, Scilab 2.6, and Tela 1.32. This time we will talk about matrices, have a look at some of the predefined functions, learn how to write our own functions, and introduce flow control statements. The article closes with a brief discussion of the applications' input and output facilities.

Matrices

Vectors help a lot if data depend on a single parameter. The different parameter values are reflected by different index values. If data depend on two parameters, vectors are a clumsy container and a more general structure, which allows for two independent indices is needed. Such a structure is called a matrix. Matrices are packed like a fresh six-pack: they are rectangular storage containers and no bottle -- oops -- element is missing.

Matrices are, for example, built from scalars as the next transcript of a GNU/Octave session demonstrates.

    octave:1> #               temperature    rain    sunshine
    octave:1> #                  degF       inches     hours
    octave:1> weather_data = [    73.4,       0.0,     10.8;  ...
    >                             70.7,       0.0,      8.5;  ...
    >                             65.2,       1.3,      0.7;  ...
    >                             68.2,       0.2,      4.1]
    weather_data =

      73.40000   0.00000  10.80000
      70.70000   0.00000   8.50000
      65.20000   1.30000   0.70000
      68.20000   0.20000   4.10000

Three new ideas appear in the example. First, we have introduced some comments to label the columns of our matrix. A comment starts with a pound sign ``#'' and extends until the end of the line. Second, the rows of a matrix are separated by semi-colons ``;'', and third, if an expression stretches across two or more lines, the unfinished lines must end with the line-continuation operator ``...''.

Similarly to vectors, matrices can not only be constructed from scalars, but from vectors or other matrices. If we had some variables holding the weather data of each day, like

    weather_mon = [73.4, 0.0, 10.8]
    weather_tue = [70.7, 0.0,  8.5]
    weather_wed = [65.2, 1.3,  0.7]
    weather_thu = [68.2, 0.2,  4.1]

we would have defined weather_data with

    weather_data = [weather_mon; weather_tue; weather_wed; weather_thu]

or, on the other hand, if we had the data from the various instruments as

    temperature = [73.4; 70.7; 65.2; 68.2]
    rain = [0.0; 0.0; 1.3; 0.2]
    sunshine = [10.8; 8.5; 0.7; 4.1]

we would have defined weather_data with

    weather_data = [temperature, rain, sunshine]

The fundamental rule is: Commas separate columns, semi-colons separate rows.

The scalars living in matrix m are accessed by applying two indices: m(row, column), where row is the row-index, and column is the column index. Thus, the amount of rain fallen on Wednesday is fetched with the expression

    octave:10> weather_data(3, 2)
    ans = 1.3000

Entries are changed by assigning to them:

    octave:11> weather_data(3, 2) = 1.1
    weather_data =

      73.40000   0.00000  10.80000
      70.70000   0.00000   8.50000
      65.20000   1.10000   0.70000
      68.20000   0.20000   4.10000

Now that we have defined weather_data we want to work with it. We can apply all binary operations that we have seen in last month's article on vectors. However, for this particular example, computing

    rain_forest_weather_data = weather_data + 2.1
    siberian_summer_weather_data = weather_data / 3.8

does not make much sense, though the computer will not complain at all. In the first example it would dutifully add 2.1 to every element of weather_data, in the second it would -- obedient like a sheepdog -- divide each element by 3.8.

Say we want to do something meaningful to weather_data and convert all temperatures from degrees Fahrenheit to degrees Celsius. To that end, we need to access all elements in the first column. The vector of interest is

    octave:16>     temp = [weather_data(1, 1); ...
    >                      weather_data(2, 1); ...
    >                      weather_data(3, 1); ...
    >                      weather_data(4, 1)]
    temp =

      73.400
      70.700
      65.200
      68.200

Obviously, the row-indices [1, 2, 3, 4] form a vector themselves. We can use a shortcut and write this vector of indices as

    temp = weather_data([1, 2, 3, 4], 1)

In general, any vector may be used as index vector. Just watch out that no index is out of range. Ordering of the indices does matter (for example weather_data([2, 1, 4, 3], 1) puts Tuesday's temperature in front) and repeated indices are permitted (for example weather_data([3, 3, 3, 3, 3, 3, 3], 1) holds Wednesday's temperature seven times).

In our example, the index-vector can be generated by a special built-in, the range generation operator ``:''. To make a vector that starts at low and contains all integers from low to high, we say

    low:high

For example

    octave:1> -5:2
    ans =

      -5  -4  -3  -2  -1   0   1   2

Our weather data example now simplifies to

    temp = weather_data(1:4, 1)

Accessing a complete column or row is so common that further shortcuts exist. If we drop both, low and high from the colon-operator, it will generate all valid column indices for us. Therefore, we reach at the shortest form to get all elements in the first column.

    octave:17> temp = weather_data(:, 1)
    temp =

      73.400
      70.700
      65.200
      68.200

With our new knowledge, we extract the sunshine hours on Tuesday, Wednesday, and Thursday

    octave:19> sunnyhours = weather_data(2:4, 3)
    sunnyhours =

      8.50000
      0.70000
      4.10000

and Tuesday's weather record

    octave:20> tue_all = weather_data(2, :)
    tue_all =

      70.70000   0.00000   8.50000

Now it is trivial to convert the data on the rain from inches to millimeters: Multiply the second column of weather_data by 25.4 (Millimeters per Inch) to get the amount of rain in metric units:

    octave:21> rain_in_mm = 25.4 * weather_data(:, 2)
    rain_in_mm =

       0.00000
       0.00000
      27.94000
       5.08000

We have already seen that vectors are compatible with scalars

    1.25 + [0.5, 0.75, 1.0]

or

    [-4.49, -4.32, 1.76] * 2

Scalars are also compatible with matrices.

    octave:1> 1.25 + [ 0.5,   0.75, 1.0; ...
    >                 -0.75,  0.5,  1.25; ...
    >                 -1.0,  -1.25, 0.5]
    ans =

      1.75000  2.00000  2.25000
      0.50000  1.75000  2.50000
      0.25000  0.00000  1.75000

    octave:2> [-4.49, -4.32, 1.76; ...
    >           9.17,  6.35, 3.27] * 2
    ans =

       -8.9800   -8.6400    3.5200
       18.3400   12.7000    6.5400

In each case the result is the scalar applied to every element in the vector or matrix.

How about vectors and matrices? Obviously, an expressions like

    [7, 4, 9] + [3, 2, 7, 6, 6]
    [2, 4; 1, 6] - [1, 1, 9, 4]

do not make any sense. In the first line the vectors disagree in size (3 vs. 5 elements), in the second line they have different shapes (2 columns and 2 rows vs. 4 columns and 1 row). To make sense, vectors or matrices that are used in an addition or subtraction must have the same shape, which means the same number of rows and the same number of columns. The technical term for ``shape'' in this context is dimension. We can query the dimension of anything with the built-in function  size().

    octave:22> size(weather_data)
    ans =

      4  3

    octave:23> size(sunnyhours)
    ans =

      3  1

The answer is a vector whose first element is the number of rows, and whose second element is the number of columns of the argument.

Multiplication and division of matrices can be defined in two flavors, both of which are implemented in the numerical workbenches.

Details

Differences

Built-In Matrix Functions

Ugh -- far too many to mention! The workbenches supply dozens of predefined functions. Here I can only wet the reader's appetite.

Generating Special Matrices
Several matrices occur often enough in computations that they have been given their own generating functions. These are for example, m-times-n matrices filled with zeros: zeros(m, n) or ones: ones(m, n), or n-times-n diagonal matrices, where the diagonal consists entirely of ones: eye(n) or the diagonal is set to numbers supplied in a vector: diag([a1, a2, ..., I<an>]).
Analyzing Matrices
Getting the smallest or largest element in matrix a: min(a), max(a), or totaling matrix a: sum(a).

Differences: GNU/Octave's min(a), max(a), and sum(a) return the column-wise result as a row vector. To get the minimum, maximum, and sum of all elements in matrix a, use min(min(a)), max(max(a)), sum(sum(a)).

Linear Algebra
We mentioned that systems of linear equations, like x * a = b, are solved for x with the slash operator ``/''. But many more linear algebra functions exist, for example singular value decomposition: svd(a), or eigenvalue computation: eig(a).

Differences: In Tela uses SVD(a) instead of svd(a), and instead of eig(a), Scilab uses spec(a) to compute the eigenvalue spectrum.

One note on performance: basically, all three applications are interpreters. This means that each expression is first parsed, then the interpreter performs desired computations, finally calling the functions inside of the expressions -- all in all a relatively slow process in comparison to a compiled program. However, functions like those shown above are used in their compiled form! They execute almost at top speed. What the interpreter does in these cases is to hand over the complete matrix to a compiled Fortran, C, or C++ function, let it do all the work, and then pick up the result.

Thus we deduce one of the fundamental rules for successful work with numerical workbenches: prefer compiled functions over interpreted code. It makes a tremendous difference in execution speed.

User Defined Functions

No matter how many functions a program may provide its users, they are never enough. Users always need specialized functions to deal with their problems, or they simply want to group repeated, yet predefined operations. In other words, there always is a need for user-defined functions.

User functions are best defined in files, so that they can be used again in later sessions. For GNU/Octave, functions files end in .m, and are loaded either automagically or with source("filename.m"). Scilab calls its function files .sci, and requires them to be loaded with getf("filename.sci"). Tela functions are stored in .t-files and loaded with source("filename.t"). As big as the differences are in loading functions, all workbenches use quite similar syntax for the definition of functions.

GNU/Octave and Scilab

    function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN)
        # function body
    endfunction

Tela

    function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN)
    {
        // function body
    };

where arg1 to argN are the functions' arguments (also known as parameters), and res1 to resN are the return values. Yes, trust your eyes, multiple return values are permitted, what might come as a surprise to most readers who are acquainted with popular programming languages. However, this is a necessity, as no function is allowed to change any of its input arguments.

Enough theory! let us write a function that takes a matrix as input and returns a matrix of the same dimensions, with the entries rescaled to lie in the interval (0, 1).

    ### Octave

    function y = normalize(x)
        ## Return matrix X rescaled to the interval (0, 1).

        minval = min(min(x))
        maxval = max(max(x))

        y = (x - minval) / (maxval - minval)
    endfunction

Now define a Scilab function that returns the spectral radius on a matrix. We use abs() which returns the magnitude of its (possibly complex) argument.

    // Scilab

    function r = spectral_radius(m)
        // Return the spectral radius R of matrix M.

        r = max(abs(spec(m)))
    endfunction

Finally, we write a Tela function which computes the Frobenius norm of a matrix.

    // Tela

    function x = frobenius(m)
    // Return the Frobenius norm X of matrix M.
    {
        x = sqrt(sum(abs(m)^2))
    };

Details:

GNU/Octave's ``automagical'' function file loading works the following way: if Octave runs into an undefined function name it searches the list of directories specified by the built-in variable LOADPATH for files ending in .m that have the same base name as the undefined function; for example, x = my_square_root(2.0) looks for the file my_square_root.m in the directories listed in LOADPATH.

Flow Control Statements

All code we have written thus far executes strictly top-to-bottom, we have not used any flow control statements such as conditionals or loops.

Before we manipulate the flow of control, we should look at logical expressions because the conditions used in conditionals and loops depend on them. Logical expressions are formed from (1.) numbers, (2.) comparisons, and (3.) logical expressions catenated with logical operators.

  1. Zero means logically false, any number not equal to zero means logically true, hence C-programmers should feel at home.
  2. The usual gang of comparison operators exist: less-than ``<'', less-or-equal ``<='', greater-than ``>'', greater-or-equal ``>='', and equal ``==''.

    Differences: The inequality operator varies quite a bit among the programs. (Octave cannot decide whether it feels like C, Smalltalk, or Pascal. Scilab wants to be Smalltalk and Pascal at the same time. :-)

        !=   ~=   <>    # Octave 
    ~=   <>         // Scilab
    !=              // Tela
  3. Complex logical expressions are formed with logical operators ``and'', ``or'' and ``not'' whose syntax is borrowed from C. However, each program uses its own set of operators. Thus, we have to list some

    Differences:

        and      or     not
    ----    ----    ----
     &       |      !  ~     # Octave
     &       |      ~        // Scilab
     &&      ||     !        // Tela

We are all set now for the first conditional, the if-statement. Note that the parenthesis around the conditions are mandatory (as they are in C). The else-branches are optional in any case.

    # Octave                // Scilab               // Tela

    if (cond)               if cond then            if (cond) {
        # then-body             // then-body            // then-body
    else                    else                    } else {
        # else-body             // else-body            // else-body
    endif                   end                     };

cond is a logical expression as described above.

while-statements:

    # Octave                // Scilab               // Tela

    while (cond)            while cond              while (cond) {
        # body                  // body                 // body
    endwhile                end                     };

Again, cond is a logical expression.

for-statements in Octave and Scilab walk through the columns of expr one by one. Most often expr will be a vector generated with the range operator ``:'', like for i = 1:10. Tela's for-statement is the same as C's.

    # Octave                // Scilab               // Tela

    for var = expr          for var = expr          for (init; cond; step) {
        # body              // body                     // body
    endfor                  end                     };

Here come some examples which only show things we have discussed so far.

Octave

    function n = catch22(x0)
        ## The famous catch-22 function: it is
        ## impossible to compute that it will
        ## stop for a specific input.  Returns 
        ## the number of loops.

        n = 0
        x = x0
        while (x != 1)
            if (x - floor(x/2)*2 == 0) 
                x = x / 2
            else
                x = 3*x + 1
            endif
            n = n + 1
        endwhile
    endfunction

Scilab

    function m = vandermonde(v)
        // Return the Vandermonde matrix M based on
        // vector V.

        [rows, cols] = size(v)
        m = []                      // empty matrix
        if rows < cols then
            for i = 0 : (cols-1)
                m = [m; v^i]
            end
        else
            for i = 0 : (rows-1)
                m = [m, v^i]
            end
        end
    endfunction

Tela

    function vp = sieve(n)
    // Sieve of Erathostenes; returns vector of
    // all primes VP that are strictly less than
    // 2*N.  1 is not considered to be a prime
    // number in sieve().
    {
        vp = #();               // empty vector
        if (n <= 2) { return };

        vp = #(2);
        flags = ones(1, n + 1);
        for (i = 0; i <= n - 2; i = i + 1)
        {
            if (flags[i + 1])
            {
                p = i + i + 3;
                vp = #(vp, p);
                for (j = p + i; j <= n; j = j + p)
                {
                    flags[j + 1] = 0
                }
            }
        }
    };

Input/Output

We have been using with the workbenches a lot. At some point we would like to call it a day, but we do not want to lose all of our work. Our functions are already stored in files. It is time to see how to make our data persist.

Simple Input and Output

All three applications at least have one input/output (I/O) model that borrows heavily from the C programming language. This model allows close control of the items read or written. Often though, it is unnecessary to take direct control over the file format written. If variables must be saved just to be restored later, simplified I/O commands will do.

Matrix Oriented I/O

As we use matrices so often, specialized functions exist to load and save whole matrices. Especially loading a matrix with a single command is convenient and efficient to read data from experiments or other programs.

Let us assume, we have the ASCII file datafile.ascii which contains the lines

    # run 271
    # 2000-4-27
    #
    # P/bar   T/K     R/Ohm
    # ======  ======  ======
    19.6      0.118352  0.893906e4
    15.9846   0.1  0.253311e5
    39.66     0.378377  0.678877e4
    13.6      0.752707  0.00622945e4
    12.4877   0.126462  0.61755e5

and sits in the current working directory. The file's five leading lines are non-numeric. They are skipped by the workbenches, but possibly aid the user in identifying her data. I have intentionally taken a data set which is not neatly formatted, as are most data files. Matrix-loading functions split the input at whitespace not at a specific column, thus they are happy with datafile.ascii.

We load the data into GNU/Octave with

    octave:1> data = load("datafile.ascii")
    data =

       1.9600e+01   1.1835e-01   8.9391e+03
       1.5985e+01   1.0000e-01   2.5331e+04
       3.9660e+01   3.7838e-01   6.7888e+03
       1.3600e+01   7.5271e-01   6.2294e+01
       1.2488e+01   1.2646e-01   6.1755e+04

or into Scilab

    -->data = fscanfMat("datafile.ascii")
     data  =

    !   19.6       0.118352    8939.06 !
    !   15.9846    0.1         25331.1 !
    !   39.66      0.378377    6788.77 !
    !   13.6       0.752707    62.2945 !
    !   12.4877    0.126462    61755.  !

or into Tela

    >data1 = import1("datafile.ascii")
    >data1
    #(      19.6,  0.118352,   8939.06;
         15.9846,       0.1,   25331.1;
           39.66,  0.378377,   6788.77;
            13.6,  0.752707,   62.2945;
         12.4877,  0.126462,     61755)

In all three examples data will contain a 5-times-3 matrix with all the values from datafile.ascii.

The complementary commands for saving a single matrix in ASCII format are

    save("data.ascii", "data")                # GNU/Octave
    fprintfMat("data.ascii", data, "%12.6g")  // Scilab
    export_ASCII("data.ascii", data)          // Tela

Note that Scilab's fprintfMat() requires a third parameter that defines the output format with a C-style template string.

Of course none of the above save commands writes the original header, the lines starting with hash-symbols, of datafile.ascii. To write these, we need the ``low-level'', C-like input/output functions, which featured in each of the three workbenches.

C-like Input/Output

For a precise control of the input and the output, C-like I/O models are offered. All three applications implement function

    printf(format, ...)

Moreover, GNU/Octave and Tela follow the C naming scheme with their C-style file I/O:

    handle = fopen(filename)
    fprintf(handle, format, ...)
    fclose(handle)

whereas Scilab prefixes these functions with an ``m'' instead of an ``f''

    handle = mopen(filename)
    mprintf(handle, format, ...)
    mclose(handle)

Whether the function is called fprintf() or mprintf(), they work the same way.

Next Month: Graphics, function plotting and data plotting.

Christoph Spiel

Chris runs an Open Source Software consulting company in Upper Bavaria/Germany. Despite being trained as a physicist -- he holds a PhD in physics from Munich University of Technology -- his main interests revolve around numerics, heterogenous programming environments, and software engineering. He can be reached at cspiel@hammersmith-consulting.com.

Copyright © 2001, Christoph Spiel.
Copying license http://www.linuxgazette.com/copying.html
Published in Issue 70 of Linux Gazette, September 2001

[ Prev ][ Table of Contents ][ Front Page ][ Talkback ][ FAQ ][ Next ]