View source for Octave Programming Tutorial/Loops and conditions

Loops are used to repeat a block of code for a known or unknown number of times, depending on the type of loop. Using loops, you will draw some nice pictures of fractals and shapes drawn with random dots.

== The <code>for</code> loop ==
We use <code>for</code> loops to repeat a block of code for a list of known values. As an example, we'll calculate the mean of a list of values. The mean is calculated from
<center><math>\overline{x} = \frac{1}{n}\sum_{i=1}^n x_i.</math></center>

We set up a vector with some values
 octave:1> x = [1.2, 6.3, 7.8, 3.6];
and calculate the mean with
 octave:2> sum = 0;
 octave:3> for entry = x,
 octave:4>   sum = sum + entry;
 octave:5> end;
 octave:6> x_mean = sum / length(x)

Line 2: Set sum equal to 0.<br />
Line 3: For each value in x, assign it to entry.<br />
Line 4: Increment sum by entry.<br />
Line 5: Ends the for loop when there are no more members of x.<br />
Line 6: Assign the final value of sum divided by the length of x to x_mean.<br />

TO DO: get a better example and explain the code.


In general, we write a <code>for</code> loop as
 for variable = vector
    ...
 end

The <code>...</code> represents the block of code that is executed exactly once for each value inside the <code>vector</code>.

=== Example: The Sierpinski triangle ===
The Sierpinski triangle is a fractal that can be generated with a very simple algorithm.
# Start on a vertex of an equilateral triangle.
# Select a vertex of the triangle at random.
# Move to the point halfway between where you are now and the selected vertex.
# Repeat from step 2.

Plotting the points that you visit by following this procedure, generates the following picture.

<center>[[Image:sierpinski_m100000.png|thumb|center|500px]]</center>

The code that generates this fractal is shown bellow. Note that this code uses one very simple for loop to generate the fractal (for i = 1:N ; ... ; end)

 axis ([-1, 1, -0.75, 1.25], "square");
 figure(1, "visible", "off");                          % no plotting window
 hold on;
 
 % defining the vertices of an equilateral triangle (symmetric to y-axis)
 V = [ 0, 1;                                           % top vertex
      cos( (pi/2)+(2*pi/3) ), sin( (pi/2)+(2*pi/3) );  % left vertex
      cos( (pi/2)-(2*pi/3) ), sin( (pi/2)-(2*pi/3) )   % right vertex
      ];
 
 r = floor(3*rand(1)+0.9999999999);                    % integer random number in [1:3]
 x = [ V(r,1), V(r,2) ];                               % initializing x on random vertex
 
 for i = 1:1000 % !!! 100000 => time=130m57.343s
   r = floor(3*rand(1)+0.9999999999);                  % integer random number in [1:3]
   x = ( x+V(r,[1:2]) )./2;                            % halfway, current to random vertex
   plot(x(1),x(2), ".");
 endfor
 
 print -dpng "sierpinski_m.png";

==== A Word of Caution ====
For performance reasons, it's best to perform as few tasks as possible inside 'for' loops, and graphical operations like plot should be moved outside of loops whenever possible. By simply storing all x values in a matrix and then plotting as shown below, the run time for the code to produce this figure drops to a few seconds, even when iterating over 100,000 elements (which the code above warns not to do).<syntaxhighlight lang="octave">
axis ([-1, 1, -0.75, 1.25], "square");
figure(1, "visible", "off");                          % no plotting window
hold on;

% defining the vertices of an equilateral triangle (symmetric to y-axis)
V = [ 0, 1;                                           % top vertex
     cos( (pi/2)+(2*pi/3) ), sin( (pi/2)+(2*pi/3) );  % left vertex
     cos( (pi/2)-(2*pi/3) ), sin( (pi/2)-(2*pi/3) )   % right vertex
     ];

r = floor(3*rand(1)+0.9999999999);                    % integer random number in [1:3]
x(1,:) = [ V(r,1), V(r,2) ];                               % initializing x as a matrix this time.

for i = 2:100000 % Safe now: 100000 => time=7.85346s
  r = floor(3*rand(1)+0.9999999999);                  % integer random number in [1:3]
  x(i,:) = ( x(i-1,:)+V(r,[1:2]) )./2;                % Add each newly calculated x value to the matrix
endfor

plot(x(:,1),x(:,2), ".");                             % plot the entire matrix just once
print -dpng "sierpinski_m.png";
</syntaxhighlight>By initializing x as a matrix of the full final size ahead of time, this processing time can be reduced still further to only 1 or 2 seconds on modern hardware. In general, if a loop can be replaced with vectorization, it should be. 

=== Exercises ===
#Write a script that sums the first <i>N</i> integers. You can check your result with the formula <math>\frac{1}{2}N(N+1)</math>.
#Write a script that does the same thing as the <code>linspace</code> function. It should start at some value, <code>xstart</code>, stop at <code>xstop</code> and create a vector that contains <i>N</i> values evenly spaced from <code>xstart</code> to <code>xstop</code>. You can use the <code>zeros</code> function to create a zero-filled vector of the right size. Use <code>help zeros</code> to find out how the function works.

== The <code>while</code> loop ==

The while loop also executes a block of code more than once but stops based on a logical condition. For example
 x = 1.0;
 while x < 1000
    disp(x);
    x = x*2;
 endwhile
will multiply <code>x</code> by 2 until its value exceeds 1000. Here, <code>x < 1000</code> is the condition of the loop. As long as the condition holds (is true), the loop will continue executing. As soon as it is false, the loop terminates and the first instruction after the loop is executed.

The general form of a while loop is

 while condition
    ...
 endwhile

=== Exercise ===
#Write a script that calculates the smallest positive integer, <i>n</i>, such that <math>a^n \ge b</math> for some real numbers <i>a</i> and <i>b</i>. (Meaning, find the smallest power of <i>a</i> that is at least <i>b</i>.) Using the <code>log</code> function is considered cheating.

=== Example: The Mandelbrot fractal ===

The Mandelbrot set is another fractal and is generated by checking how long it takes a complex number to become large. For each complex number, <i>c</i>,

#Start with <math>z_0 = 0</math>.
#Let <math>z_i = z_{i-1}^2 + c \quad\forall i = 1, 2, \ldots</math>
#Find the first <i>i</i> such that <math>|z_i| > 2</math>.

We record all of these <i>i</i> values and assign a colour to each of them. This is used to generate an image like this one.

<center>[[Image:Mandelbrot.png]]</center>

You can download the code that generates this fractal from [[Mandelbrot.m]]. Note that there is a while loop (inside some for loops) that tests whether the complex number <i>z</i> has modulus less than 2:

 while (count < maxcount) & (abs(z) < 2)
    ...
 endwhile

The first condition in the while loop checks that we do not perform too many iterations. For some values of <i>c</i> the iteration will go on forever if we let it.

See also another version by [http://nullprogram.com/blog/2007/09/02/ Christopher Wellons]

==The <code>do...until</code> statement==

These loops are very similar to while loops in that they keep executing based on whether a given condition is true or false. There are however some important difference between <code>while</code> and <code>do...until</code> loops.
# <code>while</code> loops have their conditions at the beginning of the loop;<br><code>do...until</code> loops have theirs at the end.
# <code>while</code> loops repeat as long as the condition is true;<br><code>do...until</code> loops continue as long as theirs is false.
# <code>while</code> will execute 0 or more times (because the condition is at the beginning);<br><code>do...until</code> loops will execute 1 or more times (since the condition is at the end).

The general form of a <code>do...until</code> loop is
 do
    ...
 until condition

=== Exercise ===
Write a script that calculates the greatest common divisor (GCD) of two positive integers. You can do this using Euclid's algorithm.

=== Challenge ===
Write a script that generates random number pairs (a, b) that are distributed uniformly
#over the disc <math>\left\{(x, y) | x^2 + y^2 \le 1\right\}</math> (the first image below);
#as in the second image below

<center>[[Image:Octave_uniform_random_circle.png]] [[Image:Octave uniform random challenge.png]]</center>

== The <code>break</code> and <code>continue</code> statements ==

Sometimes it is necessary to stop a loop somewhere in the middle of its execution or to move on to the next value in a for loop without executing the rest of the loop code for the current value. This is where the <code>break</code> and <code>continue</code> statements are useful.

The following code demonstrates how the break statement works.

 total = 0;
 while true
    x = input('Value to add (enter 0 to stop): ');
    if x == 0
       break;
    endif
    total = total+x;
    disp(['Total: ', num2str(total)]);
 endwhile

Without the <code>break</code> statement, the loop would keep executing forever since the condition of the <code>while</code> loop is always true. The <code>break</code> allows you to jump past the end of the loop (to the statement after the endwhile).

The <code>break</code> statement can be used in any loop: <code>for</code>, <code>while</code> or <code>do...until</code>.

The continue statement also jumps from the inside of a loop but returns to the beginning of the loop rather than going to the end. In a
# <code>for</code> loop, the next value inside the vector will be assigned to the for variable (if there are any left) and the loop restarted with that value;
# <code>while</code> loop, the condition at the beginning of the loop will be retested and the loop continued if it is still true;
# <code>do...until</code> loop, the condition at the end of the loop will be tested and the loop continued from the beginning if it is still false.

As an example, the following code will fill the lower triangular part of a square matrix with 1s and the rest with 0s.

 N = 5;
 A = zeros(N); % Create an N x N matrix filled with 0s
 
 for row = 1:N
    for column = 1:N
       if column > row
          continue;
       endif
       A(row, column) = 1;
    endfor
 endfor
 
 disp(A);

Note that the inner <code>for</code> skips (continues) over the code that assigns a 1 to an entry of <code>A</code> whenever the column index is greater than the row index.

== The <code>if</code> statement ==

The general form of the <code>if</code> statement is

 if condition1
    ...
 elseif condition2
    ...
 else
    ...
 endif

If <code>condition1</code> evaluates to true, the statements in the block immediately following the <code>if</code> are executed. If <code>condition1</code> is false, the next condition (<code>condition2</code> in the <code>elseif</code>) is checked and its statements executed if it is true. You can have as many <code>elseif</code> statements as you like. The final set of statements, after the <code>else</code>, is executed if all of the conditions evaluate to false. Note that the <code>elseif</code> and <code>else</code> parts of the <code>if</code> statement are optional.

The following are all valid <code>if</code> statements:

 % Take the log of the absolute value of x
 if x > 0
    y = log(x);
 elseif x < 0
    y = log(-x);
 else
    disp("Cannot take the log of zero.");
 endif
 
 x = input("Enter a value: ");
 if x > 0
    disp("The number is positive");
 endif
 
 if x < 0
    disp("The number is negative");
 endif
 
 if x == 0
    disp("The number is zero");
 endif

=== Example: The fractal fern ===

This algorithm is not quite complete. Have a look at the .m file available from http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4372&objectType=file.

The image to the right can be generated with the following algorithm:

   1. Let x1 and y1 be random values between 0 and 1.
   2. Choose one of the linear transformations below to calculate (xi+1, yi+1) from (xi, yi):
         1. xi+1 = 0
            yi+1 = 0.16yi
         2. xi+1 = 0.20xi − 0.26yi
            yi+1 = 0.23xi + 0.22yi + 1.6
         3. xi+1 = −0.15xi + 0.28yi
            yi+1 = 0.26xi + 0.24yi + 0.44
         4. xi+1 = 0.85xi + 0.04yi
            yi+1 = −0.04xi + 0.85yi + 1.6
      The first transformation is chosen if probability 0.01, the second and third with probability 0.07 each and the fourth with probability 0.85.
   3. Calculate these values for i up to at least 10,000.

You can download the code that generates this fractal as fracfern.m (this is disabled for now).

----
Return to the [[Octave Programming Tutorial|Octave Programming tutorial index]]

{{BookCat}}