Thursday, 5 September 2013

Metaprogramming in R with an example: Beating lazy evaluation

Functional languages allows us to treat functions as types. This brings us a distinct advantage of being able to write a code that generates further code, this practise is generally known as metaprogramming. As a functional language R project provides tools to perform well structured code generation. In this post, I will present a simple example that generates functions on the fly based on different parametrisation in the function body. Consider the following simple function taking a vector as an argument and returning the number of element that are higher than a given threshold. 
1
2
3
4
myFun <- function(vec) {
  numElements <- length(which(vec > threshold))
  numElements
}
If somehow we need to have a different threshold value within the body, for a moment accept that it is a requirement rather than proposing to have an other argument in the function definition. Instead of rewriting the function by hand we will write a function that generates all these functions in our work space. Problematic bit of this exercise will be to beat lazy evalution.  Here is the function that produces losts of myFun type functions:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
genMyFuns <- function(thresholds) {
  ll <- length(thresholds)
  print("Generating functions:")
  for(i in 1:ll) {
    fName <- paste("myFun.", i, sep="")
    print(fName)
    assign(fName, eval(
                       substitute(
                                  function(vec) {
                                    numElements <- length(which(vec > tt));
                                    numElements;
                                  }, 
                                  list(tt=thresholds[i])
                                 )
                      ),
             envir=parent.frame()
           )
  }
}
Let's shortly analyse  this function. If we don't use substitute explicitly there, due to lazy evalution our value for the threshold will not be assigned at the loop value but the  last value of thresholds[i]. Here is one numeric example on the R CLI session:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>  genMyFuns(c(7, 9, 10))
[1] "Generating functions:"
[1] "myFun.1"
[1] "myFun.2"
[1] "myFun.3"
>  myFun.1(1:20)
[1] 13
>  myFun.2(1:20)
[1] 11
>  myFun.3(1:20)
[1] 10
> 
To be able to generate code is very powerful tool. However, a caution should be taken in practicing code generation in a large project. This may bring more problems in debugging. Every powerful method comes with a hidden cost.

Thursday, 13 June 2013

Practicing static typing in R: Prime directive on trusting our functions with object oriented programming

The creator of S language which R is derived from John Chambers said in one of his books  Software for data analysis programming with R
...This places an obligation on all creators of software to program in such a
way that the computations can be understood and trusted. This obligation
I label the Prime Directive.
He was referring to prime directive from Star Trek. One of the practice in this direction is to have a proper checks in place for the types we use. We can trust that if we pass for example a wrong type to our function, it will fail gracefully. So a type system of a programming language is quite important in mission critical numerical computations. Since R language is weakly typed language, or dynamically typed similar to Perl, Python or Matlab/Octave, most of R users omit to place type checks in their functions if not rarely. For example take the following function that takes arguments of a matrix, a vector and a function name. It applies the named function to each columns of the matrix listed in the given vector. Assuming named function is returning a single number our function will return a vector of numbers.

myMatrixOperation  <-  function(A, v, fName) {
  sliceA <-  A[, v];
  apply(sliceA, 2, fName);
}

One obvious way to put if statements for each argument in our function. So, function may look like:

myMatrixOperation <- function(A, v, fName) {
  if(!is.matrix(A)) {
   stop("A is not a matrix");
  }
  if(!is.vector(v)) {
   stop("v is not a vector");
  }
  if(!is.funcion(fName)) {
   stop("fName is not a function");
  }
  sliceA <- A[, v];
  apply(sliceA, 2, fName);
}
The problem with this approach appears to be the fact that it is too verbose and if we have a repeating pattern of arguments in many functions and many arguments, we would copy and paste code many times. It would not only look ugly but wastes our time. Luckily there is a mechanism to address this: S4-class system. Let's define an S4 class for our set of arguments, following an example instantiation.

setClass("mySlice", representation(A="matrix", v="vector", fName="function"))

myS <- new("mySlice",, A=matrix(rnorm(9),3,3),v=c(1,2), fName=mean)
str(myS)
Formal class 'mySlice' [package ".GlobalEnv"] with 3 slots
  ..@ A    : num [1:3, 1:3] 0.356 -0.34 -0.642 -0.466 2.915 ...
  ..@ v    : num [1:2] 1 2
  ..@ fName:function (x, ...)

Now if we re-write the function that uses our S4 class with type checking only to passing object once.
is.mySlice <- function(obj) {
  l <- FALSE 
  if(class(obj)[1] == "mySlice") { l <- TRUE } 
  l 
} 

myMatrixOperation <- function(mySliceObject) { 
  if(!is.mySlice(mySliceObject)) { 
    stop("argument is not class of mySlice") 
  }  
  sliceA <- mySliceObject@A[, mySliceObject@v]; 
  apply(sliceA, 2, mySliceObject@fName); 
} 
This simple example demonstrates how we can introduce a good organization to our R codes, that obeys the prime directive. Further more modern approach to object orientation is introduced by John Chambers called Reference classes. If you practice this kind of approach in your R codes than I can only say; Live long and prosper.

Monday, 6 May 2013

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)

Wednesday, 3 April 2013

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.
(c) Copyright 2008-2024 Mehmet Suzen (suzen at acm dot org)

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License