File: spiel.html

package info (click to toggle)
lg-issue70 1-3
links: PTS
area: main
in suites: woody
size: 1,228 kB
ctags: 146
sloc: makefile: 36; perl: 29; sh: 4
file content (1167 lines) | stat: -rw-r--r-- 36,314 bytes
<!--startcut  ==============================================-->
<!-- *** BEGIN HTML header *** -->
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML><HEAD>
<title>Numerical Workbenches, part II LG #70</title>
</HEAD>
<BODY BGCOLOR="#FFFFFF" TEXT="#000000" LINK="#0000FF" VLINK="#0000AF"
ALINK="#FF0000">
<!-- *** END HTML header *** -->

<CENTER>
<A HREF="http://www.linuxgazette.com/">
<IMG ALT="LINUX GAZETTE" SRC="../gx/lglogo.png" 
	WIDTH="600" HEIGHT="124" border="0"></A> 
<BR>

<!-- *** BEGIN navbar *** -->
<IMG ALT="" SRC="../gx/navbar/left.jpg" WIDTH="14" HEIGHT="45" BORDER="0" ALIGN="bottom"><A HREF="mcgucken.html"><IMG ALT="[ Prev ]" SRC="../gx/navbar/prev.jpg" WIDTH="16" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="lg_toc70.html"><IMG ALT="[ Table of Contents ]" SRC="../gx/navbar/toc.jpg" WIDTH="220" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><A HREF="../lg_frontpage.html"><IMG ALT="[ Front Page ]" SRC="../gx/navbar/frontpage.jpg" WIDTH="137" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="http://www.linuxgazette.com/cgi-bin/talkback/all.py?site=LG&article=http://www.linuxgazette.com/issue70/spiel.html"><IMG ALT="[ Talkback ]" SRC="../gx/navbar/talkback.jpg" WIDTH="121" HEIGHT="45" BORDER="0" ALIGN="bottom"  ></A><A HREF="../lg_faq.html"><IMG ALT="[ FAQ ]" SRC="./../gx/navbar/faq.jpg"WIDTH="62" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="tranter.html"><IMG ALT="[ Next ]" SRC="../gx/navbar/next.jpg" WIDTH="15" HEIGHT="45" BORDER="0" ALIGN="bottom"  ></A><IMG ALT="" SRC="../gx/navbar/right.jpg" WIDTH="15" HEIGHT="45" ALIGN="bottom">
<!-- *** END navbar *** -->
<P>
</CENTER>

<!--endcut ============================================================-->

<H4 ALIGN="center">
"Linux Gazette...<I>making Linux just a little more fun!</I>"
</H4>

<P> <HR> <P> 
<!--===================================================================-->

<center>
<H1><font color="maroon">Numerical Workbenches, part II</font></H1>
<H4>By <a href="mailto:cspiel@hammersmith-consulting.com">Christoph Spiel</a></H4>
</center>
<P> <HR> <P>  

<!-- END header -->




<p>In 
<A HREF="../issue69/spiel.html">Part&nbsp;I</A> ,
we looked at the most basic operations of the numerical
workbenches GNU/Octave&nbsp;2.1.34, Scilab&nbsp;2.6, and Tela&nbsp;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.</p>

<h2><a name="matrices">Matrices</a></h2>

<p>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.</p>

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

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

<pre>
      73.40000   0.00000  10.80000
      70.70000   0.00000   8.50000
      65.20000   1.30000   0.70000
      68.20000   0.20000   4.10000
</pre>

<p>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&nbsp;``<code>#</code>'' and extends until the end of the line. Second,
the rows of a matrix are separated by semi-colons&nbsp;``<code>;</code>'', and
third, if an expression stretches across two or more lines, the unfinished
lines must end with the line-continuation operator ``<code>...</code>''.</p>

<p>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</p>

<pre>
    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]
</pre>

<p>we would have defined <code>weather_data</code> with</p>

<pre>
    weather_data = [weather_mon; weather_tue; weather_wed; weather_thu]
</pre>

<p>or, on the other hand, if we had the data from the various instruments
as</p>

<pre>
    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]
</pre>

<p>we would have defined <code>weather_data</code> with</p>

<pre>
    weather_data = [temperature, rain, sunshine]
</pre>

<p>The fundamental rule is: <em>Commas separate columns, semi-colons separate
rows.</em></p>

<p>The scalars living in matrix&nbsp;<code>m</code> are accessed by applying
two indices: <code>m(row, column)</code>, where <em>row</em> is the row-index,
and <em>column</em> is the column index. Thus, the amount of rain fallen on
Wednesday is fetched with the expression</p>

<pre>
    octave:10&gt; weather_data(3, 2)
    ans = 1.3000
</pre>

<p>Entries are changed by assigning to them:</p>

<pre>
    octave:11&gt; weather_data(3, 2) = 1.1
    weather_data =
</pre>

<pre>
      73.40000   0.00000  10.80000
      70.70000   0.00000   8.50000
      65.20000   1.10000   0.70000
      68.20000   0.20000   4.10000
</pre>

<p>Now that we have defined <code>weather_data</code> 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</p>

<pre>
    rain_forest_weather_data = weather_data + 2.1
    siberian_summer_weather_data = weather_data / 3.8
</pre>

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

<p>Say we want to do something meaningful to <code>weather_data</code> 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</p>

<pre>
    octave:16&gt;     temp = [weather_data(1, 1); ...
    &gt;                      weather_data(2, 1); ...
    &gt;                      weather_data(3, 1); ...
    &gt;                      weather_data(4, 1)]
    temp =
</pre>

<pre>
      73.400
      70.700
      65.200
      68.200
</pre>

<p>Obviously, the row-indices&nbsp;<code>[1,&nbsp;2,&nbsp;3,&nbsp;4]</code>
form a vector themselves. We can use a shortcut and write this vector of
indices as</p>

<pre>
    temp = weather_data([1, 2, 3, 4], 1)
</pre>

<p>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 <code>
weather_data([2, 1, 4, 3], 1)</code> puts Tuesday's temperature in front) and
repeated indices are permitted (for example <code>weather_data([3, 3, 3, 3, 3,
3, 3], 1)</code> holds Wednesday's temperature seven times).</p>

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

<pre>
    low:high
</pre>

<p>For example</p>

<pre>
    octave:1&gt; -5:2
    ans =
</pre>

<pre>
      -5  -4  -3  -2  -1   0   1   2
</pre>

<p>Our weather data example now simplifies to</p>

<pre>
    temp = weather_data(1:4, 1)
</pre>

<p>Accessing a complete column or row is so common that further shortcuts
exist. If we drop both, <em>low</em> and <em>high</em> 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.</p>

<pre>
    octave:17&gt; temp = weather_data(:, 1)
    temp =
</pre>

<pre>
      73.400
      70.700
      65.200
      68.200
</pre>

<p>With our new knowledge, we extract the sunshine hours on Tuesday,
Wednesday, and Thursday</p>

<pre>
    octave:19&gt; sunnyhours = weather_data(2:4, 3)
    sunnyhours =
</pre>

<pre>
      8.50000
      0.70000
      4.10000
</pre>

<p>and Tuesday's weather record</p>

<pre>
    octave:20&gt; tue_all = weather_data(2, :)
    tue_all =
</pre>

<pre>
      70.70000   0.00000   8.50000
</pre>

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

<pre>
    octave:21&gt; rain_in_mm = 25.4 * weather_data(:, 2)
    rain_in_mm =
</pre>

<pre>
       0.00000
       0.00000
      27.94000
       5.08000
</pre>

<p>We have already seen that vectors are compatible with scalars</p>

<pre>
    1.25 + [0.5, 0.75, 1.0]
</pre>

<p>or</p>

<pre>
    [-4.49, -4.32, 1.76] * 2
</pre>

<p>Scalars are also compatible with matrices.</p>

<pre>
    octave:1&gt; 1.25 + [ 0.5,   0.75, 1.0; ...
    &gt;                 -0.75,  0.5,  1.25; ...
    &gt;                 -1.0,  -1.25, 0.5]
    ans =
</pre>

<pre>
      1.75000  2.00000  2.25000
      0.50000  1.75000  2.50000
      0.25000  0.00000  1.75000
</pre>

<pre>
    octave:2&gt; [-4.49, -4.32, 1.76; ...
    &gt;           9.17,  6.35, 3.27] * 2
    ans =
</pre>

<pre>
       -8.9800   -8.6400    3.5200
       18.3400   12.7000    6.5400
</pre>

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

<p>How about vectors and matrices? Obviously, an expressions like</p>

<pre>
    [7, 4, 9] + [3, 2, 7, 6, 6]
    [2, 4; 1, 6] - [1, 1, 9, 4]
</pre>

<p>do not make any sense. In the first line the vectors disagree in size (3
vs. 5&nbsp;elements), in the second line they have different shapes
(2&nbsp;columns and 2&nbsp;rows vs. 4&nbsp;columns and 1&nbsp;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&nbsp;<code>
size()</code>.</p>

<pre>
    octave:22&gt; size(weather_data)
    ans =
</pre>

<pre>
      4  3
</pre>

<pre>
    octave:23&gt; size(sunnyhours)
    ans =
</pre>

<pre>
      3  1
</pre>

<p>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.</p>

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

<ul>
<li>Element by element multiplication or division of two vectors or matrices
of same dimensions: The number in the first row and first column of the first
matrix is multiplied by the number in the first row and first column of the
second matrix and so on for every element. 

<pre>
    a = [3, 3; ...
         6, 4; ...
         6, 3]
    b = [9, 3; ...
         8, 2; ...
         0, 3]
</pre>

<pre>
    octave:1&gt; a .* b
    ans =
</pre>

<pre>
      27   9
      48   8
       0   9
</pre>

<p>The element-by-element operators are preceded by a dot: element-by-element
multiplication ``<code>.*</code>'' and element-by-element division
``<code>./</code>''.</p>
</li>

<li>Matrix multiplication as known from Linear Algebra: <em>
c</em>&nbsp;=&nbsp;<em>a</em>&nbsp;* <em>b</em>, where <em>a</em> is a <em>
p</em>-times-<em>q</em>&nbsp;matrix and <em>b</em> is a <em>
q</em>-times-<em>r</em>&nbsp;matrix. The result&nbsp;<em>c</em> is a <em>
p</em>-times-<em>r</em> matrix. 

<p>Example:</p>

<pre>
    a = [3, 3; ...
         6, 4; ...
         6, 3]
</pre>

<pre>
    b = [-4,  0, 1, -4; ...
         -1, -3, 2,  0]
</pre>

<pre>
    octave:1&gt; a * b
    ans =
</pre>

<pre>
      -15   -9    9  -12
      -28  -12   14  -24
      -27   -9   12  -24
</pre>

<p>Although we have not seen <code>for</code>-loops yet (they will be discussed
<a href="#flow control statements">farther down</a>), I would like to write
the code behind the matrix multiplication operator&nbsp;``<code>*</code>'' to
give the reader an impression of the operations involved.</p>

<pre>
    for i = 1:p
        for j = 1:r
            sum = 0
            for k = 1:q
                sum = sum + a(i, k)*b(k, j)
            end
            c(i, j) = sum
        end
    end
</pre>

<p>Compare these triply nested <code>for</code>-loops with the simple
expression <code>c&nbsp;=&nbsp;a&nbsp;*&nbsp;b</code>.</p>
</li>

<li>Matrix division? You cannot divide by a matrix! However, operator
``<code>/</code>'' is defined for vectors and matrices. But writing <em>
x</em>&nbsp;= <em>b</em>&nbsp;/&nbsp;<em>a</em>, where <em>a</em> and <em>
b</em> are matrices or vectors has nothing to do with division at all! It
means: please solve the system of linear equations 

<pre>
    x * a = b
</pre>

<p>for <em>x</em>, given matrix&nbsp;<em>a</em> and the
right-hand-side(s)&nbsp;<em>b</em>. Here ``<code>*</code>'' denotes matrix
multiplication as defined in the previous item, and the same rules for
compatible dimensions of <em>a</em> and <em>b</em> apply.</p>

<pre>
    a = [-2, 3,  1; ...
          7, 8,  6; ...
          2, 0, -1]
</pre>

<pre>
    b = [-26,  5, -6; ...
          24, 53, 26]
</pre>

<pre>
    octave:1&gt; x = b / a
    x =
</pre>

<pre>
       7.00000  -2.00000   1.00000
       7.00000   4.00000   5.00000
</pre>

<p>Isn't that an easy way to solve a system of linear equations? Imagine you
had to write the code which does exactly that.</p>

<p>Finally, let us verify the result by multiplying with <em>a</em> again</p>

<pre>
    octave:2&gt; x*a
    ans =
</pre>

<pre>
      -26.0000    5.0000   -6.0000
       24.0000   53.0000   26.0000
</pre>

<p>which, as expected, recovers <em>b</em>.</p>
</li>
</ul>

<p><strong>Details</strong></p>

<ul>
<li>For convenience GNU/Octave and Scilab define an alternative matrix
division operator&nbsp;``<code>\</code>''. <em>x</em>&nbsp;=&nbsp;<em>
a</em>&nbsp;\&nbsp;<em>b</em> solves the linear system of equations 

<pre>
    a * x = b
</pre>

<p>for <em>x</em>, given matrix&nbsp;<em>a</em> and the
right-hand-side(s)&nbsp;<em>b</em>. This is the form most users prefer,
because here <em>x</em> is a column vector, whereas
operator&nbsp;``<code>/</code>'' returns <em>x</em> as row-vector.</p>
</li>

<li>operator&nbsp;``<code>\</code>'' has the dotted cousin ``<code>.\</code>''
and the relation <em>a</em>&nbsp;./&nbsp;<em>b</em>&nbsp;==&nbsp;<em>
b</em>&nbsp;.\&nbsp;<em>a</em> holds.</li>
</ul>

<p><strong>Differences</strong></p>

<ul>
<li>Scilab and Tela use C++-like comments 

<pre>
    // This is a Scilab or a Tela comment
</pre>
</li>

<li>Tela does not need or understand the line continuation operator
``<code>...</code>'' 

<pre>
    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)
</pre>

<p>In interactive mode, Tela does not handle multi-line expressions as the
above. Multi-line expressions must be read from a file (with <code>
source("filename.t")</code>).</p>
</li>

<li>In Tela the operators ``<code>*</code>'' and ``<code>/</code>'' work
element by element, this is, they work like ``<code>.*</code>'' and
``<code>./</code>'' do in GNU/Octave and Scilab. Matrix multiplication
(<em>a</em>&nbsp;*&nbsp;<em>b</em> in GNU/Octave or Scilab) is written as 

<pre>
    a ** b
</pre>

<p>or</p>

<pre>
    matmul(a, b)
</pre>

<p>solving systems of linear equations (<em>b</em>&nbsp;/&nbsp;<em>a</em> in
Octave or Scilab) as</p>

<pre>
    linsolve(a, b)
</pre>
</li>
</ul>

<h2><a name="builtin matrix functions">Built-In Matrix Functions</a></h2>

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

<dl>
<dt><strong><a name="item_Generating_Special_Matrices">Generating Special
Matrices</a></strong><br>
</dt>

<dd>Several matrices occur often enough in computations that they have been
given their own generating functions. These are for example, <em>
m</em>-times-<em>n</em> matrices filled with zeros: <code>zeros(m, n)</code>
or ones: <code>ones(m, n)</code>, or <em>n</em>-times-<em>n</em> diagonal
matrices, where the diagonal consists entirely of ones: <code>eye(n)</code> or
the diagonal is set to numbers supplied in a vector: <code>diag([a1, a2, ...,
I&lt;an&gt;])</code>.</dd>

<dt><strong><a name="item_Analyzing_Matrices">Analyzing
Matrices</a></strong><br>
</dt>

<dd>Getting the smallest or largest element in matrix&nbsp;<em>a</em>: <code>
min(a)</code>, <code>max(a)</code>, or totaling matrix&nbsp;<em>a</em>: <code>
sum(a)</code>. 

<p><strong>Differences:</strong> GNU/Octave's <code>min(a)</code>, <code>
max(a)</code>, and <code>sum(a)</code> return the column-wise result as a row
vector. To get the minimum, maximum, and sum of all elements in
matrix&nbsp;<em>a</em>, use <code>min(min(a))</code>, <code>
max(max(a))</code>, <code>sum(sum(a))</code>.</p>
</dd>

<dt><strong><a name="item_Linear_Algebra">Linear Algebra</a></strong><br>
</dt>

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

<p><strong>Differences:</strong> In Tela uses <code>SVD(a)</code> instead of
<code>svd(a)</code>, and instead of <code>eig(a)</code>, Scilab uses <code>
spec(a)</code> to compute the eigenvalue spectrum.</p>
</dd>
</dl>

<p>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.</p>

<p>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.</p>

<h2><a name="user defined functions">User Defined Functions</a></h2>

<p>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.</p>

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

<p>GNU/Octave and Scilab</p>

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

<p>Tela</p>

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

<p>where <em>arg1</em> to <em>argN</em> are the functions' arguments (also
known as parameters), and <em>res1</em> to <em>resN</em> 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.</p>

<p>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&nbsp;(0,&nbsp;1).</p>

<pre>
    ### Octave
</pre>

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

<pre>
        minval = min(min(x))
        maxval = max(max(x))
</pre>

<pre>
        y = (x - minval) / (maxval - minval)
    endfunction
</pre>

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

<pre>
    // Scilab
</pre>

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

<pre>
        r = max(abs(spec(m)))
    endfunction
</pre>

<p>Finally, we write a Tela function which computes the Frobenius norm of a
matrix.</p>

<pre>
    // Tela
</pre>

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

<p><strong>Details:</strong></p>

<a name="automagical_explanation"></a>
<p>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 <code>LOADPATH</code> for files
ending in .m that have the same base name as the undefined function; for
example, <code>x = my_square_root(2.0)</code> looks for the file <em>
my_square_root.m</em> in the directories listed in <code>LOADPATH</code>.</p>

<h2><a name="flow control statements">Flow Control Statements</a></h2>

<p>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.</p>

<p>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.)&nbsp;numbers,
(2.)&nbsp;comparisons, and (3.)&nbsp;logical expressions catenated with
logical operators.</p>

<ol>
<li>Zero means logically false, any number not equal to zero means logically
true, hence C-programmers should feel at home.</li>

<li>The usual gang of comparison operators exist: less-than
``<code>&lt;</code>'', less-or-equal ``<code>&lt;=</code>'', greater-than
``<code>&gt;</code>'', greater-or-equal ``<code>&gt;=</code>'', and equal
``<code>==</code>''. 

<p><strong>Differences:</strong> 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. :-)</p>

<pre>
    !=   ~=   &lt;&gt;    # Octave 
    ~=   &lt;&gt;         // Scilab
    !=              // Tela
</pre>
</li>

<li>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 

<p><strong>Differences:</strong></p>

<pre>
    and      or     not
    ----    ----    ----
     &amp;       |      !  ~     # Octave
     &amp;       |      ~        // Scilab
     &amp;&amp;      ||     !        // Tela
</pre>
</li>
</ol>

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

<pre>
    # Octave                // Scilab               // Tela
</pre>

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

<p><em>cond</em> is a logical expression as described above.</p>

<p><code>while</code>-statements:</p>

<pre>
    # Octave                // Scilab               // Tela
</pre>

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

<p>Again, <em>cond</em> is a logical expression.</p>

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

<pre>
    # Octave                // Scilab               // Tela
</pre>

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

<p>Here come some examples which only show things we have discussed so
far.</p>

<p>Octave</p>

<pre>
    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.
</pre>

<pre>
        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
</pre>

<p>Scilab</p>

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

<pre>
        [rows, cols] = size(v)
        m = []                      // empty matrix
        if rows &lt; 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
</pre>

<p>Tela</p>

<pre>
    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 &lt;= 2) { return };
</pre>

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

<h2><a name="input/output">Input/Output</a></h2>

<p>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.</p>

<h3><a name="simple input and output">Simple Input and Output</a></h3>

<p>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.</p>

<ul>
<li>Octave offers the most flexible solution with the <code>
save</code>/<code>load</code> command pair. 

<pre>
    save filename varname1 varname2 ... varnameN
</pre>

<p>saves the variables named <em>varname1</em>, <em>varname2</em>, ..., <em>
varnameN</em> in file&nbsp;<em>filename</em>. The complementary</p>

<pre>
    load filename varname1 varname2 ... varnameN
</pre>

<p>command restores them from <em>filename</em>. If <code>load</code> is given
no variable names, all variables form <em>filename</em> are loaded. Handing
over names to <code>load</code> selects only the named variables for
loading.</p>

<p>Note that the <code>save</code> and <code>load</code> commands do not have
parenthesis and their arguments are separated by spaces not commas. Filename
and variable names are strings.</p>

<pre>
    save "model.oct-data" "prantl" "reynolds" "grashoff"
    load "model.oct-data" "reynolds"
</pre>

<p>By default <code>load</code> does not overwrite existing variables, but
complain with an error if the user tries to do so. When it is save to discard
of the values of existing variables, add option&nbsp;``<code>-force</code>''
to <code>load</code>, like</p>

<pre>
    load -force "model.oct-data" "reynolds"
</pre>

<p>and variable&nbsp;<code>reynolds</code> will be loaded from file <em>
model.oct-data</em> no matter whether it has existed before or not.</p>
</li>

<li>Scilab's simple I/O parallels that of GNU/Octave: 

<pre>
    save(filename, var1, var2, ..., varN)
</pre>

<p>However, the variables <em>var1</em>, ..., <em>varN</em> are not strings,
but appear literally. This means that the name of a variable is not stored in
the file. The association between name and contents is lost!</p>

<p>The complementary function</p>

<pre>
    load(filename, varname1, varname2, ..., varnameN)
</pre>

<p>restores the contents of <em>filename</em> in the variables named <em>
varname1</em>, <em>varname2</em>, ... <em>varnameN</em>.</p>
</li>

<li>Tela lets the users save her variables with the 

<pre>
    save(filename, varname1, varname2, ..., varnameN)
</pre>

<p>function, preserving the association between variable name and variable
contents. The complementary</p>

<pre>
    load(filename)
</pre>

<p>function loads all variables stored in <em>filename</em>. It is not
possible to select specific variables.</p>
</li>
</ul>

<h3><a name="matrix oriented i/o">Matrix Oriented I/O</a></h3>

<p>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.</p>

<p>Let us assume, we have the ASCII file <em>datafile.ascii</em> which
contains the lines</p>

<pre>
    # 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
</pre>

<p>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
<em>datafile.ascii</em>.</p>

<p>We load the data into GNU/Octave with</p>

<pre>
    octave:1&gt; data = load("datafile.ascii")
    data =
</pre>

<pre>
       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
</pre>

<p>or into Scilab</p>

<pre>
    --&gt;data = fscanfMat("datafile.ascii")
     data  =
</pre>

<pre>
    !   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.  !
</pre>

<p>or into Tela</p>

<pre>
    &gt;data1 = import1("datafile.ascii")
    &gt;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)
</pre>

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

<p>The complementary commands for saving a single matrix in ASCII format
are</p>

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

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

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

<h3><a name="clike input/output">C-like Input/Output</a></h3>

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

<pre>
    printf(format, ...)
</pre>

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

<pre>
    handle = fopen(filename)
    fprintf(handle, format, ...)
    fclose(handle)
</pre>

<p>whereas Scilab prefixes these functions with an ``<code>m</code>'' instead
of an ``<code>f</code>''</p>

<pre>
    handle = mopen(filename)
    mprintf(handle, format, ...)
    mclose(handle)
</pre>

<p>Whether the function is called <code>fprintf()</code> or <code>
mprintf()</code>, they work the same way.</p>

<p><EM>Next Month: Graphics, function plotting and data plotting.</EM></p>





<!-- *** BEGIN bio *** -->
<SPACER TYPE="vertical" SIZE="30">
<P> 
<H4><IMG ALIGN=BOTTOM ALT="" SRC="../gx/note.gif">Christoph Spiel</H4>
<EM>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 
<A HREF="mailto:cspiel@hammersmith-consulting.com">cspiel@hammersmith-consulting.com</A>.</EM>

<!-- *** END bio *** -->

<!-- *** BEGIN copyright *** -->
<P> <hr> <!-- P --> 
<H5 ALIGN=center>

Copyright &copy; 2001, Christoph Spiel.<BR>
Copying license <A HREF="../copying.html">http://www.linuxgazette.com/copying.html</A><BR> 
Published in Issue 70 of <i>Linux Gazette</i>, September 2001</H5>
<!-- *** END copyright *** -->

<!--startcut ==========================================================-->
<HR><P>
<CENTER>
<!-- *** BEGIN navbar *** -->
<IMG ALT="" SRC="../gx/navbar/left.jpg" WIDTH="14" HEIGHT="45" BORDER="0" ALIGN="bottom"><A HREF="mcgucken.html"><IMG ALT="[ Prev ]" SRC="../gx/navbar/prev.jpg" WIDTH="16" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="lg_toc70.html"><IMG ALT="[ Table of Contents ]" SRC="../gx/navbar/toc.jpg" WIDTH="220" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><A HREF="../lg_frontpage.html"><IMG ALT="[ Front Page ]" SRC="../gx/navbar/frontpage.jpg" WIDTH="137" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="http://www.linuxgazette.com/cgi-bin/talkback/all.py?site=LG&article=http://www.linuxgazette.com/issue70/spiel.html"><IMG ALT="[ Talkback ]" SRC="../gx/navbar/talkback.jpg" WIDTH="121" HEIGHT="45" BORDER="0" ALIGN="bottom"  ></A><A HREF="../lg_faq.html"><IMG ALT="[ FAQ ]" SRC="./../gx/navbar/faq.jpg"WIDTH="62" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="tranter.html"><IMG ALT="[ Next ]" SRC="../gx/navbar/next.jpg" WIDTH="15" HEIGHT="45" BORDER="0" ALIGN="bottom"  ></A><IMG ALT="" SRC="../gx/navbar/right.jpg" WIDTH="15" HEIGHT="45" ALIGN="bottom">
<!-- *** END navbar *** -->
</CENTER>
</BODY></HTML>
<!--endcut ============================================================-->