QuickTip: Adding Search Paths in Octave

Suppose you want a utility (such as matrixTeX.m) to be visible to every octave session.

1. Select a suitable folder to keep your local octave files in (e.g. /home/sachin/octave/ on my Linux box). If such a folder doesn't exist, make one.

2.  Fire up an octave session

3. In Octave, add the path to the path variable:

octave:> addpath(genpath("/home/sachin/octave")) 

4. List directories in default path (your output will likely look different)

octave:> path 

/home/sachin/octave
/home/sachin/octave/io-2.4.11
/home/sachin/octave/io-2.4.11/doc
/home/sachin/octave/io-2.4.11/packinfo
/home/sachin/octave/io-2.4.11/templates

etc.

4. So far, we haven't permanently written the modified path to disk. To do this:

octave:> savepath

Floating Point Numbers are Not Associative

On a recent exam, I asked students for suggestions on how to test for equality of two floating point numbers.

Loren Shure and Mike Croucher present simple examples which highlight the perils. For example, (0.1 + 0.2) + 0.3 is not equal to 0.1 + (0.2 + 0.3).

As pointed out in a comment, this also implies:

octave:2> x = 0.3 + 0.2 + 0.1;
octave:3> y = 0.1 + 0.2 + 0.3;
octave:4> x==y
ans = 0
octave:5> x-y
ans = -1.1102e-16

Nicer Octave Plots for Presentations/Documents

GNU Octave has a very capable plotting system (which is transitioning away from gnuplot as its backend, even as we "speak"). Quite often, you generate a figure and make it look like you want:

You then want to incorporate it into some other document - say a presentation or a document.

So you say:

print -dpng test0.png

Note that you could also use -depsc2 or -dpdf to generate EPS or PDF formats among a whole host of other possibilities.

You get a figure that looks like:

Clearly imperfect.

There are a number of ways to improve the image.

These include writing to tikz, TeX, which gives you exquisite control over all features, but this can take too much effort.

Or you could spend time fiddling with adjustable parameters to make the plot look just right. If you are trying to work on "the" central plot of a paper, or a PhD thesis, then such care is certainly warranted.

But if you are merely dissatisfied (as I am) with the size and font of the text, there is an easy fix.

print -dpng -FHelvetica:18 test1.png

Conditional Column Switching in Octave or Matlab

Consider a simple Matlab/Octave exercise.

I have an N by 2 vector x. Thus, x(i,1) and x(i,2) are the entries of the i-th row. Say I want to rearrange x so that the first entry in any row is larger than the second. That is, I want to switch the entries of a row, only if a certain condition is met.

It is easy to write this in a element-by-element fashion

for i = 1:N
  if x(i,2) > x(i,1)
    c = x(i,1);
    x(i,1) = x(i,2);
    x(i,2) = c;
  endif
endfor

But we all know that element-by-element operations in Matlab or Octave suck! Here is another vectorized method:


SwitchCond = x(:,2) < x(:,1);
x(SwitchCond,[1,2]) = x(SwitchCond,[2,1]);

Besides being faster, it is also more succinct.

Making "journal quality" graphs in Matlab

A very nice blog on making journal quality graphs by extensively using handles at Loren and the Art of Matlab.

As a side note, the actual blog entry was written by guest blogger Jiro Doke, who I went to school with, as mentioned earlier.

Setting up the random number generator seed in Matlab and Octave

Mike Croucher as an interesting post on correctly and incorrectly setting up the seed in Matlab.

The key takeaway is that one should not use rand('seed',x) or rand('state',x) to reset the random number seed, since Matlab defaults to an inferior generator used in older versions. It is preferable to use rng(x) instead.

GNU Octave on the other hand has a less quirky take on the whole thing.

Using rand('seed',x) defaults to older generators, while rand('state',x) resets the seed, but still uses a version of Mersenne Twister.

Cleve Moler: Matlab, textbooks and blog

As many of you who use Matlab might know, Cleve Moler was a professor of math and computer science at a bunch of universities (including Michigan), before he co-founded Mathworks.

I should have known this before, given how much Matlab and GNU Octave I use, but only recently did I find out about two free textbooks that he wrote. They are available (as pdfs - either chapterwise or in full) on the Mathworks website. The first called "Numerical Computing with MATLAB" is an engaging trek through the usual topics of numerical analysis with Matlab. It is written in an extremely accessible style, and highlights the many special insights he has accumulated over the years. The use of Matlab allows him to quickly explore some topics, do some quick computations, plot and visualize, and mine for knowledge and wisdom.

The other book called "Experiments with Matlab" is similar in spirit, except that it is targeted at younger audiences. Rather than using a standard numerical analysis course to lay down the skeleton of the book, it jumps from one puzzle to another, one curiosity to the next looking at the underlying patterns and math.

Both books are very highly recommended, regardless of whether you are a student learning this material for the first time, or an instructor who has taught similar topics several times. There is something in it for everybody.

I should also mention that Cleve writes an engaging blog called "Cleve's Corner" which I now have on my Google Reader.

Non-Uniform Quadrature Points in Octave/Matlab

Suppose you are given a bunch of points in the form of a vector x (between "a" and "b") and an accompanying density function p(x)>0. For concreteness, consider

x = [0:0.1:1.0]'; px = 0.3 - (x-0.5).^2;

The top subfigure depicts the density function p(x), and the cumulative density function C(x). The bottom subfigure depicts the  quadrature points "z" (circles), and the lines mark the end of the intervals, using N=11.

This density function is shown in the picture above.

We want to create a N*1 vector "z" of quadrature points that are distributed according to the density p(x). We also want a vector  "dz" corresponding to the width of strip belonging to a particular quadrature point. As an additional constraint, we assume that we want z(1) = a, and z(N)  = b to match the end-points.

The program attached at the end is able to do this relatively efficiently. In fact, I used the following command to generate the picture above:

[z dz] = GridDensity(x,px,11,1);

GNU Octave Code:

%

%  PROGRAM:

% Takes in a PDF or density function, and spits out a bunch of points in

%       accordance with the PDF

%

%  INPUT:

%       x  = vector of points. It need *not* be equispaced

%       px = vector of same size as x: probability distribution or

%            density function. It need not be normalized but has to be positive.

%   N  = Number of points >= 3. The end points of "x" are included

%       necessarily

%       Pt = Optional argument. If present, then some plotting.

% 

%  OUTPUT:

%       z  = Points distributed according to the density

%       hz = width of the "intervals" - useful to apportion domain to points

%            if you are doing quadrature with the results, for example.

%

%  (*) Sachin Shanbhag, March 5, 2012

%  GNU Public License

%

function [z h] = GridDensity(x,px,N,Pt)

  npts = 100;                              % can potentially change

  xi   = linspace(min(x),max(x),npts)';    % reinterpolate on equi-spaced axis

  pint = interp1(x,px,xi,'spline');        % smoothen using splines

  ci   = cumtrapz(xi,pint);                

  pint = pint/ci(npts);

  ci   = ci/ci(npts);                      % normalize ci

  alfa = 1/(N-1);                          % alfa/2 + (N-1)*alfa + alfa/2

  zij  = zeros(N,1);                       % quadrature interval end marker

  z    = zeros(N,1);                       % quadrature point

  z(1)  = min(x);  

  z(N)  = max(x); 

%

% ci(Z_j,j+1) = (j - 0.5) * alfa

%

  beta  = [0.5:1:N-1.5]'*alfa;

  zij   = [z(1); interp1(ci,xi,beta,'spline'); z(N)];

  h        = diff(zij);

  clear beta;

%

% Quadrature points are not the centroids, but rather the center of masses

% of the quadrature intervals

%

  beta     = [1:1:N-2]'*alfa;

  z(2:N-1) = interp1(ci,xi,beta,'spline');

%

% Some plotting if required

%

  if(nargin>3)

    subplot(2,1,1)

    plot(xi,ci,'b-','LineWidth',2,xi,pint,'r-','LineWidth',2);

    axis([min(xi),max(xi)]);

    legend('CDF','PDF')

    title('Visualization')

    xlabel('x');

    ylabel('CDF(x)/PDF(x)')

    subplot(2,1,2)

    plot(z,0.5*ones(size(z)),'ro');

    axis([min(xi),max(xi)]);

    xlabel('z');

    hold on;

    for i = 2:N

      X = [zij(i); zij(i)];

      Y = [0; 1];

      plot(X,Y,'b')

    endfor

    hold off;

  endif

endfunction

Nested Subfunctions in Matlab and Octave

In a previous post, I mentioned subfunctions in Octave and Matlab.

Nested functions are an interesting new beast. I can see situations where they can be useful as explained in this short video, and this longer, more detailed blog.

They are currently not featured in GNU Octave. Here's their reasoning:

There was discussion in Octave of having these even prior to Matlab, and the decision was made not to have these in Octave at the time for compatibility.

Now that Matlab has recently introduced nested functions, Octave will probably have them soon as well. Until then nested functions in Octave are treated as sub-functions with the same scoping rules as sub-function.

The authors of Octave consider the nested function scoping rules of Matlab to be more problems than they are worth as they introduce diffiult to find bugs as inadvertantly modifying a variable in a nested function that is also used in the parent is particularly easy.

Writing multiple column vectors to file in Matlab and Octave

Consider a simple example. If x and y are two vectors of the same size, one obvious way of writing them to file is:

n  = 1e5;
x  = 1:n;
y  = x.^2;
f1 = fopen('xy.dat','w');
for i = 1:length(x)
    fprintf(f1,'%d\t%d\n',x(i),y(i));
end
fclose(f1);

[Side Note: Matlab does not have a command like Octave's "printf" - if "fprintf" is used without a file handle, then it outputs to the screen.]

I ran timings for the "red" writing part of the code above in Matlab 7.4 R2011b and Octave 3.4.2 on my Desktop, and found that it took 3.4s and 2.6s, respectively. It isn't usual for Octave to outrun Matlab, but this may have something to do with the peculiarities of my installation.

I replaced the "for-loop" in the write operation with a single command:

fprintf(f1,'%d\t%d\n',x,y);

and my timings improved. It took 0.13s and 0.11s for Matlab and Octave, respectively. So clearly getting rid of the for loop helps (no surprise!).

Sometimes, you may be in a hurry and not care so much for the control over formatting that fprintf offers. If you want avoid all the fopen and fclose business, you could instead use, something like:

xy = [x' y'];
save('xy.dat', 'xy','-ascii');

Note that I need to define a "matrix", with the columns side by side to use this.  It took 0.23s and 0.28s for Matlab and Octave, respectively - so there is a small performance hit.