Octave/Matlab: Loop through Files in a Directory

Looping through files that match a certain pattern is an easy and routine operation for a BASH shell for example:

One can mimic this feature from within an Octave or Matlab program by using the command "glob". Here's how:

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

Passing parameters in Matlab or Octave: fzero and quad

Consider a simple function f(x) = 2x written in Matlab/Octave (myfunc.m).

function y = myfunc(x)

    y = 2*x

end

You could easily use the functions fzero and quad, to find a root of f(x) and the definite integral of f(x) between limits, say "a" and "b".

You would say something like:

fzero('myfunc',1.0)

where 1.0 is the guessed root. The answer that Matlab finds is indeed zero.

or,

quad('myfunc',0.0,1.0)

where a=0.0, and b=1.0 are the lower and upper limits. The quad function computes the correct answer (1.0).

Let us say that we had a slightly more general function f(x) = c*x, where "c" is a parameter (c=2, in the previous case). That is, we have the Matlab/Octave function:

function y = myfunc(x,c)

    y = c*x

end

Quite often we want to use fzero and quad on such a function for a particular setting of "c". That is we would like to be able to pass the parameter "c" to the function.

There are several ways to doing this. Some of them (like making "c" a global variable are ugly). The one I find most convenient is the following:

c=2;

quad(@(x) myfunc(x,c), 0, 1.0)

or,

quad(@(x) myfunc(x,2.0), 0,1.0)

fzero works in similar fashion. You could say:

fzero(@(x) myfunc(x,2.0), 1.0)

Just to reiterate, this method works in both Matlab and Octave, and you could easily generalize it to passing multiple parameters.

Profiling Fortran and C programs

Given how easy it is to profile programs written in C and Fortran, I wonder why many people don't do it. Profiling lets you monitor the performance of your code.

One way to do it is to use GPROF.

Here is how you use it.

1. Compile your code (test.f90) with the -pg flag

gfortran -pg test.f90 (or gcc -pg test.c)

2. Run the program

a.out

This creates a report gmon.out, which is unfortunately not in readable format.

3. To read it, say

gprof a.out

You will see a reasonably detailed report on number of function calls and time spent in different parts of the program.

Orthogonal Polynomials: Mathematica

Orthogonal polynomials are everywhere.

The following Mathematica program takes in a weight function (as a function of x), the domain (a and b), and spits out the first "n" corresponding orthonormal polynomials.

OrthoPoly[w_, a_, b_, n_] := Module[{monoBasis, T},
  monoBasis = Table[{x^i}, {i, 0, n - 1}];
  oP = monoBasis;
  oP[[1]] = oP[[1]]/Sqrt[Integrate[w*oP[[1]]*oP[[1]], {x, a, b}]];
  For[i = 2, i <= n, i = i + 1,
   For[T = 0; j = 1, j < i, j = j + 1,
    T = T + Integrate[w*oP[[i]]*oP[[j]], {x, a, b}]*oP[[j]];
    ];
   oP[[i]] = oP[[i]] - T;
   oP[[i]] =
    oP[[i]]/Sqrt[Integrate[w*oP[[i]]*oP[[i]], {x, a, b}]] //
     Simplify;
   ] ;
  oP
  ]

Here's a screenshot, for the first four Chebyshev and Legendre polynomials (click to enlarge). Note that these polynomials are unique up to a multiplicative constant. Orthonormality freezes that constant.

Numerical Differentiation using Mathematica

I am teaching myself Mathematica over summer, and coded up a simple module that computes numerical differentiation formulas automatically. Essentially, the module enables one to compute arbitrary differentiation rules like those explained here.

This module spits out the different formulae for a m-point approximation to the n-th derivative, and the leading error term.

Here is the (updated) module:

AppDeriv[m_, n_] := Module[{points, IP, nthDeriv, e1},
  points = Table[{x0 + i h, Subscript[f, i]}, {i, 0, m - 1}];
  IP = InterpolatingPolynomial[points, x];
  nthDeriv = D[IP, {x, n}];
  e1 = (1/Factorial[m])*
    D[Apply[Times, (x - Table[{x0 + i h}, {i, 0, m - 1}])], {x, n}];
  Formula =
   TableForm[
    Table[{Subsuperscript[f, i, n], nthDeriv /. x -> (x0 + i h),
       Superscript[f, m] e1 /. x -> (x0 + i h)}, {i, 0, m - 1}] //
     Simplify ]]

This seems to work alright for me (click to enlarge).

Dispersing points on the surface of a sphere

Dispersing or choosing points, uniformly or randomly, on the surface of a sphere is a task that someone like me who does molecular simulations, runs into, every once in a while. In some form or shape, this problem is encountered, for example, while modeling a 3D random walk, the motion of a Brownian particle in space, decorating a nanoparticle with stabilizers etc.

While the task is not hard, there are some elegant methods and potential pitfalls that a person doing this for the first time should be aware of.

As I mentioned earlier, one may choose points randomly, or uniformly on the surface (of a unit sphere, here, but can be generalized trivially).

1. Random:

The key point here is that it is incorrect to choose points by selecting angles phi and theta corresponding to spherical coordinates from uniform distributions. This incorrectly concentrates points near the poles.

The correct method is to select two random numbers u and v from a uniform distribution, and construct the two angles as follows:

u = rand();            // uniform random number between

v = rand();            // zero and one.

phi = 2 * pi * u;      // pi is 3.141...

theta = acos(2*v - 1)  // acos is inverse cosine

This is explained with figures here.

2. Uniform:

This is trickier than I thought when I first attempted to do it. Here is an old link (circa 1998, but still good!) which touches upon some of the subtleties. The trouble starts with what "uniform" means, for an arbitrary number of points to be distributed.

If we decide that uniform means a distribution which "maximizes the minimum separation" between n points on the sphere, we can start going somewhere.

There are a number of algorithms that can try to solve this problem, as described on this excellent page. Here is C++ code that I wrote while implementing the golden spiral rule from that page.

//
// Spray n points uniformly on the surface of a sphere 

// of radius "radius"
//

// The positions are returned as an array of Vectors

// where, struct Vector { double x; double y; double z }

//

// Based on algorithm from 

// http://www.xsi-blog.com/archives/115

//

void SprayPointsSphere(int n, double radius, Vector * p)
{
  double inc = PI * (3.0 - sqrt(5.0));
  double off = 2.0/n;

  for(int k = 0; k < n; k++) {
    double y   = k * off - 1 + (off/2);
    double r   = sqrt(1 - y*y);
    double phi = k * inc;

    p[k].x = cos(phi)*r * radius;

    p[k].y = y * radius;

    p[k].z = sin(phi)*r * radius);

  }
}