Showing posts with label MATLAB. Show all posts
Showing posts with label MATLAB. Show all posts

Sunday, 8 September 2013

Weighted Histograms: MATLAB or Octave

Histograms are widely used in exploratory data analysis.  In some applications, such as Monte Carlo simulations, sometimes the data or the distribution manifest with corresponding weights, for example in the analysis of  umbrella sampling results. In those situations weights of each data point should be taken into account in constructing histograms out of data points.  If you are programming with R,  weights package would help you to do this. Recently, I have added small utility function in mathworks file exchange called Generate Weighted Histograms to plot weighted histograms. In this post, I would like to demonstrate how these function can be used. Let's generate some synthetic data, which contains values and corresponding weights and find the weighted histogram.

1
2
3
4
5
6
7
% Generate synthetic data
myData.values = rand(100,1); 
myData.weights = rand(100,1);
% Generate weighted histogram
[histw, vinterval] = histwc(myData.values, myData.weights, 10);
% Visualize
bar(vinterval, histw) 
This trivial example demonstrate the usage of histwc function. One can adjust the bins as the third argument of the function.


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)

Friday, 1 February 2013

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.

Monday, 3 December 2012

Function closures and abstraction in MATLAB/Octave

Recently a short idea on how function closures can be build in R is shown [link]. It is also discussed in detail by Darren Wilkinson's blog [link]. Here we will mimic the approach presented there. Recall that function closures are basically a reference to a function with the attachment of its lexical environment. Note that, it is not a function pointer while it also return the values of its arguments if they are defined in the given scope. This sort of language construction can be very handy when some of the inner details need to be parametrized, so to speak to prevent re-writing the function over and over again.

A trivial example of adding a number to a given number can be formulated.  Let's say x is a number and we need to add 4 to x.

function add4(x)
  x+4

end

So far so good, but do you really want to re-write an other function to add 7. Not really! So instead we can generate a function that generates addition function for 7 (Sounds like an onion.). This is quite standard in functional languages and noting new actually!

>> addx = @(x) @(y) x + y ;
>> add4=addx(4);
>> add7=addx(7);

>> add4(5)

>> add7(5)
12

So function add4 and add7 are generated by the function addx. We just utilized function handles. This might be quite trivial. Let's try the technique in a little more realistic setting. Imagine we need to take a median of each vector, which has different sizes stored in a structure.  Let's generate one

>> rng(0.42,'v4'); % reproduce numbers
>> myStr.vec1 =  rand(11,1);
>> myStr.vec2 =  rand(5,1);
>> myStr.vec3 =  rand(20,1);
>> myStr.vec4 =  rand(8,1);

So a function that takes a structure like this and returns the median value of the named member

>> medianVec= @(str, memName) median(str.(memName));

For example for vec1

>> medianVec(myStr, 'vec1');
 0.5194

If we need this procedure for vec1 repetitively, we don't need to re-write the second argument all the time. So

>> medianVecN = @(memName)  @(str) median(str.(memName));
>> medianVec1 = medianVecN('vec1');

>> medianVec1(myStr); % Vala!
0.5194 

This might still look trivial but hopefully it would give you an idea of the technique.


MATLAB/OCTAVE crash tutorial for programmers: Useful bits

In this post, I would like to give you a crash tutorial on the basics of  MATLAB ™ - GNU Octave (MO). They are not fully compatible but quite close in grammar. For differences see here. We will concentrate on the basic data structures and object manipulation that would help you to write useful code quickly.
  • Vectors/matrices and Indexing
    Vectors or vectorial operations are the core of MO
    >> A = [ 3 6 7 8 2];
    >> A(1:3); % range indexing will print values from index 1 to index 3
    >> A(4:end); % Last two values of this vector
    >> A > 4 % logical indexing; this will return logical vector, components satisfying the condition
    0 1 1 1 0
    >> A(A > 4)  % one can use this logical indexing directly to extract elements that satisfy the condition
    6 7 8

    A([2 3 4]) % positional indexing
    logInx = A>4
    A(logInx)
    6 7 8
    Let's define 3 by 3 matrix
    >> M = [ 1 3 4 ; 5 6 9 ; 2 3 4]
    >> M(:,1:2) % first two columns

  • Structures and Cells
    A structure definition can be performs as follows
    persons=struct();
    Fields can be defined after .
    persons.age=12; Further index can be assigned dynamically
    persons(2).age=32;
    persons(3).age=40; 
  • Functions
    Return vector with  function definition must be in a seperate file. For example here pov.m (in the working directory)

    function [ pov2, pov3] = pov(x)
      pov2= x^2; pov3 = x^3;
    end


    (pwd will tell you which directory you are working in. A return vector here is [pov2, pov3].

  • Function handles and apply functions
    To define new function ff with function handle
    it takes a structure xx and returns its field age value

    ff = @(xx) xx.age

    We can apply this function to each element of the structure,
    returning an array (vector)

    arrayfun(ff, persons)

  • Example Task:
    Lets align two vectors, meaning that we found matching indices and delete non matching ones. Here is the solution:

    p1 = [2    11    15    19    28    41    45    47    48    51];
    p2 = [2    11    31    36    41    45    51    60    68    71];
    [~,indexp1] = intersect(p1,p2);
    % if you use ~, it means don't assign anything
    [~,indexp2] = intersect(p2,p1);
    p1 = p1(indexp1);
    p2 = p2(indexp2);

Friday, 26 October 2012

Intersect many sets in MATLAB efficiently

Finding intersection of two sets in MATLAB is one of the most frequently used function if you are doing lots of data merging tasks.  This is achieved by intersect function. However if you would like to do the same for many sets,  there are not many alternatives.

Consider that you placed your arrays into a cell vector and would like to get all intersections of a given vector among other remaining arrays and its corresponding indices uniquely. The following function will do the job for you  efficiently. It is not O($N^2$) algorithm, so computational cost only grows linearly. I called this function intersectall, code is maintained in the mathworks file exchange; intersectall.


Wednesday, 10 October 2012

Matlab/Octave: imagesc with variable colorbar

Plotting a discrete matrix with a colour map is a common task that one would need in many different fields of sciences. To plot a discrete matrix with MATLAB, a possibility is to use imagesc, while it handles discrete points better then other density utilities, where it can handle (0,0) placement better. However, if there are limited number of values in the matrix entries, you may want to keep only existing ones on the color labels. One way to achieve this is to re-write your 'colormap' via 'jet', that maps existing colours on the current matrix. However, this is not sufficient, you need to re-write your matrix as well. See example code below. So called a variable 'colormap' situation is handled with this approach. (This is tested with MATLAB, on Octave, I think one needs to be careful on 'colormap', values behaviour of 'jet' is little different.)

% bare minimum - example -
clear all
close all
allLevelNames   = {'one', 'two', 'three', 'four'};
levelValues     = [1 2 3 4];
myMatrix        =  [3 3 3 3; 3 3 3 3; 4 4 1 1; 4 4 1 1];
% Now re-write colormap
aV              = length(levelValues);
jj              = jet(aV);
availableValues = unique(myMatrix);
getIndexes      = arrayfun(@(xx) find(levelValues == xx), availableValues); 
colormap(jj(getIndexes,:)); 
% Re-write myMatrix in the availableValues range 
myMatrix = arrayfun(@(x) find(availableValues == x), myMatrix);
% Now  Plot
h    = imagesc(myMatrix);

% Handle if there is more then one color is in use
maxV = max(availableValues); 
minV = min(availableValues); 
if(length(availableValues) > 1) caxis([minV maxV]); end; 
hcb=colorbar; 
set(hcb, 'YTick', availableValues); 
set(hcb, 'YTickLabel', allLevelNames(getIndexes));
set(gca,'YDir','normal') % Not to invert x-y but putting (0,0) on the bottom left 
xlabel('x values');
ylabel('y values (care on (0,0) location from imagesc)');
title({['Variable colormap handling'];['by msuzen at gmail']});





(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