Logarithmic Density Plots MATLAB or Octave: Two Gaussians Different Scales

Common visualization of 3D data manifest as surface and mesh plots. Data may contain values which are different in several orders of magnitude if not more. In this post, I will demonstrate how to transform your 3D data values to logarithmic scale so that this magnitude discrepancy is not washed out due to linear scaling. I will give code sample using MATLAB but example should be applicable to Octave as well.
But the procedure should be generic and easily be written in any other language of your choice.

The most common 3D data set appear as values over spatial dimensions. For example third dimension could be a probability distribution. Let's say we use two multivariate gaussians with different means and covariances using mvnpdf.  In the first figure we see only a small region; visiable for us. However there are data points which are washed out due to different scale. In the code below we re-write the data's zeros to a very small value (the matrix F), say 1-e14, so that we can convert the values to logarithmic scale,  remember that log(0) is an -Inf and not plotable. Resulting figure gives us more information about our data.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

% Two Gaussians on the plane mu1 = [0 0]; sigma1 = [.25 .3; .3 1]; mu2 = [1, 1]; sigma2 = [.0025 .003; .003 0.1]; x = -4:.1:4; y = -4:.1:4; [X, Y] = meshgrid(x,y); F = mvnpdf([X(:) Y(:)], mu1, sigma1) + mvnpdf([X(:) Y(:)], mu2, sigma2); F = reshape(F,length(x),length(y)); figure(1) pcolor(X, Y, F); shading flat xlabel('x'); ylabel('y'); c=colorbar; ylabel(c,'Probability Density') ; figure(2) % Fill zeros with small number smallN = 1e-14; [~,y] = size(F); getZeroIds = @(colN) find(F(:,colN) < smallN); applyAllCols = @(ii) evalin('base', ['colN=', sprintf('%d',ii) ,';F(getZeroIds(colN), colN)=smallN;']); arrayfun(applyAllCols, 1:y); % pcolor(X, Y, log10(F)); shading flat xlabel('x'); ylabel('y'); c=colorbar; ylabel(c,'Probability Density') ; % clear all close all [x y] = meshgrid(-5:0.5:5, -5:0.5:5); z = - (1 - x.^2 -y.^2) .* exp(-(x.^2 + y.^2)).* exp(-(0.5*x.^2 - y.^2)); surf(z);

Figure 1: Two Gaussians in linear scale.

Figure 2: Two Gaussians in logarithmic scale
(see the  code above)

Reading Binary Data Files Written in C: Python - Numpy Case Study

One of the major components of  scientific computation is producing results into a file. Usually bunch of numbers, usually structured, is written into so called data files. Writing  binary files as outputs is practised commonly. It is the choice due to good speed and/or memory efficiency compare to plain ASCII files. However, a care on documenting  endianness must be observed for good portability.

There are many popular languages to achieve the above task. For speed and efficiency reasons usually C or Fortran is used in writing out a binary file. Let's give an example of writing one integer (42) and three doubles (0.01, 1.01, 2.01) into binary file in C:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

#include <stdio.h> #include <stdlin.h> int main() { FILE * myF; int i,j; double *numbers, kk; myF = fopen("my.bindata", "wb") ; numbers = malloc(3*sizeof(double)); i = 42; fwrite(&amp;i, sizeof(int), 1, myF); for(j=0; j<3; j++) { kk = (double)j+1e-2; numbers[j] = kk; } fwrite(numbers, sizeof(double), 3, myF); fclose(myF); return(0); }

This code would produce a binary file called my.bindata. Our aim is to read this into Python so we can post-process the results i.e. visualisation or further data analysis. The core idea is to use higher language in processing the outputs directly instead of writing further C code; so to speak avoiding one more step in our work flow and avoiding cumbersome compilation of extra C code.

In order to read from files byte by byte, the standard library of Python provides a module called struct. Basically this module provides packing and unpacking of data into or from binary sources, in this case study our source is a file. However it is tedious and error prone to use this in a custom binary file where format would contain different types. Well at least needs an effort to read our custom binary file. At this point, our friend is Numpy facilities. Specially two functionality;
numpy.dtype and numpy.fromfile. The former provides an easy way of defining our file's format similar to Fortran syntax via creation of a data type object as its name stands. The later is a direct way of reading the binary file in one go that would return us a Python object that contains the all information present in the data file.
Here is the Numpy code that reads our binary file created by the above C code.

1 2 3 4 5 6 7 8 9 10 11

import numpy as np dt = np.dtype("i4, (3)f8") myArray = np.fromfile('my.bindata', dtype=dt) myArray[0] #(42, [0.01, 1.01, 2.01]) myArray[0][1] #array([ 0.01, 1.01, 2.01]) myArray[0][0] #42 myArray[0][1][1] #1.01

 
I have tested this case study on GNU/Linux PC, so the binary file is little-endian hence the writing and reading patterns.  Ideally a generic wrapper around this Python code would help to simplify things.

Big-O notation: Care on asymptotic behaviour

Big-O notation conveys us information about how a given algorithm run time scales with respect to a given input size.  However, one has to be careful on the interpretation of this scaling (limiting) behaviour. First of all, when we talk about $O(N^2)$, it doesn't mean that it is valid for any input size. Let's  illustrate that asymptotic behaviour shown by big-O notation could be “misleading” for “smaller” input. Let’s take number of operation $T(N)=4 N^2 – 2N$. $T’(4)$ must be 4 times of $T’(2)$ if we consider as $T’(N) = N^2$. But if we use $T(N)$ instead, $T(4)$ would be $4.666667$ times of $T(2)$.

A good example in this direction  is N-body force calculations (or here). An $O(N^2)$ algorithm performs better compare to $O(NlogN)$ ones up to a break-even point even 50K particles. This is a good paper explaining this;

Pringle, Gavin J. "Comparison of an O (N) and an O (N log N) N-body solver."

[pdf].

Multicore Run in Matlab via Python: Generation of Henon Maps

Hennon Map for a=1.4 and b=0.3.

Multi-core processors are considered as a de facto standard. In scientific computing it is common to face a problem of generating data based on a single functionality with many different parameters. This can be thought as "single instruction multiple parameter sets" type of computation. It is quite natural to utilise all cores available in your host machine. While many people uses MATLAB for rapid prototyping; Here I show how to generate many Hennon Maps with different initial conditions (ICs) using MATLAB and Python. We will use Python to drive the computations and "single instruction" is being a function in MATLAB.

Let's first shortly remember the definition of the Hennon Map

$ x_{n+1} =  y_{n} + 1 - \alpha x_{n}^2 $
$ y_{n+1} = \beta x_{n}  $

It is known that for parameter ranges $\alpha \in [1.16, 1.41]$ and $\beta \in [0.2, 0.3]$ map generates chaotic behaviour i.e. sensitive dependence on initial conditions.

Here is a MATLAB function that generates a png file for a given parameters, initial condition values, file name and the upper bound for the iteration.

function hennonMap(a, b, xn, yn, upper, filename)
% HENNONMAP generate hennon map at given parameters and initial conditions
%
% Mehmet Suzen
% msuzen on gmail
% (c) 2013
% GPLv3
%
% a \in [1.16 1.41]
% b \in [0.2, 0.3]
%
% Example of running this in the command line
%  > matlab  -nodesktop -nosplash -logfile my.log -r "hennonMap(1.4, 0.3, -1.0, -0.4, 1e4, 'hennonA1.4B0.3.png'); exit;"
%
  Xval = zeros(upper, 1);
  Yval = zeros(upper, 1);
  for i=1:upper
    Xval(i) = xn;
    Yval(i) = yn;
    x  = yn + 1 - a*xn^2;
    y  = b * xn;
    xn = x;
    yn = y;
  end
  h = figure
  plot(Xval, Yval,'o')
  title(['Henon Map at ', sprintf(' a=%0.2f', a), sprintf(' b=%0.2f', b)])
  xlabel('x')
  ylabel('y')
  print(h, '-dpng', filename)
end

Running this function from a command line, as described in the function help, would generate the figure shown above. Now imagine that we need to generate many plots with different ICs. It is easy to  open many shells and run from command line many times. However it is a manual task and we developers do not like that at all!  Python would help us in this case, specifically its multiprocessing module, to spawn the process from a single script. The concept here is closely related to threading, while multiprocessing module mimics threading API.

For example, let's say we have 16 cores and we would like to generate Henon maps for 16 different initial conditions. Here is a Python code running the above Hennon Map function by invoking 16 different MATLAB instances with each using different ICs:

#
# Mehmet Suzen
# msuzen on gmail
# (c) 2013
# GPLv3
#  
#  Run Hennon Map on multi-process (spawning processes)
#

from multiprocessing import Pool
import commands
import numpy
import itertools

def f(argF):
        a            = '%0.1f' % argF[0]
        b            = '%0.1f' % argF[1]
        filename     = "hennonA" + a + "B" + b + ".png"
        commandLine  = "matlab  -nodesktop -nosplash -logfile my.log -r \"hennonMap(1.4, 0.3," +  a + "," + b + ", 1e4, '" + filename + "'); exit;\""
        print commandLine
        return commands.getstatusoutput(commandLine)

if __name__ == '__main__':
        pool     = Pool(processes=12)              # start 12 worker processes
        xns      = list(numpy.linspace(-1.0, 0.5, 4))
        yns      = list(numpy.linspace(-0.4, 1.1, 4))
        print pool.map(f, list(itertools.product(xns, yns)))

It would be interesting to do the same thing only using R. Maybe in the next post, I'll show that too.