## Random number generation with different probability distribution – MATLAB

When programming with MATLAB, often you need to generate random numbers either uniformly or with a normal distribution. There are two main functions, namely rand and randn, that take care of this problem. Here I want to show you a deeper way to generate numbers, not only with the above mentioned distributions also with other distributions such as: weibul, logistic, loglogistic, gamma, nakagami, binomial, negative binomial, poisson, rician, etc.

1- using a function called makedist you can make a distribution by providing it with corresponding variables of the specific distribution. For instance the value for variance for a normal distribution.

2- using a function called random you can then generate as much as numbers you need from the distribution made in the previous step.

Here I show you some example:

pd_uniform = makedist('Uniform');
mean(pd_uniform)
std(pd_uniform)
uniform_rand_numbers = random(pd_uniform,10000,2);
subplot(1,3,1);scatter(uniform_rand_numbers(:,1),uniform_rand_numbers(:,2),'.');
title('Uniform distribution [0 1]');

pd_normal = makedist('Normal')
mean(pd_normal)
std(pd_normal)
normal_rand_numbers = random(pd_normal,10000,2);
subplot(1,3,2);scatter(normal_rand_numbers(:,1),normal_rand_numbers(:,2),'.');
title('Normal distribution \sigma = 1, \mu = 0');

pd_gamma = makedist('Gamma','a',1,'b',1)
mean(pd_gamma)
std(pd_gamma)
gamma_rand_numbers = random(pd_gamma,10000,2);
subplot(1,3,3);scatter(gamma_rand_numbers(:,1),gamma_rand_numbers(:,2),'.');
title('Gamma distribution a = 1 , b = 1');


## Check the version when a function was introduced – MATLAB

I just saw the following function written in MATLAB that can extract the information about when  a specific function was introduced by MATLAB. Check it out here.

## How to carry out an operator with probability p – Optimization Algorithms

When implementing optimization algorithms, we often bump into this phrase: “an operator needs to be carried out with probability P”. For instance, this can happen when implementing the crossover or mutation phases in Evolutionary Algorithms. It means that you want to select an operation with a specific probability (chance). Let’s have an example: If we want to increase a variable, x, by one unit with probability of p = 0.2. This means that 20% of times we want to carry out this operation and 80% of times we prefer not to. To achieve this behavior, the easiest way is to use a uniform distribution and sample a value in range [0 1]. If the sample is less than the probability p then we can carry out the operation. In MATLAB, the code would be as follows:

x = 0;
p = 0.2;
if rand < p
x = x + 1;
end


## MEX file in MATLAB – Part 1

In this tutorial I will show you a very simple example of a MEX file in MATLAB. A MEX file is a compiled version of a C/C++ code written and compiled in MATLAB that can be used/executed as an M-file. A MEX file has a way better performance respect to a m-file in MATLAB. Thus if you want to increase the performance of your code, converting it into a MEX file is one of the ways.
1 – To create your first experience open MATLAB and create a multiplier.c file in the workspace.
2- The first thing you should add to your code is header files. The most useful one in MATLAB is mex.h, so add it to your code using the include command:

#include "mex.h"


3- In this example and generally in any code you want to write you will have some input/output for the program and you have to check and pass them into the program. When using MEX files, you can write a gateway function for this purposes:

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *inVector1;              /* input matrix */
double *inVector2;               /* 1xN input matrix */
size_t ncols;                   /* size of matrix */
double *product;              /* output scalar */

/* create a pointer to the real data in the input matrix  */
inVector1 = mxGetPr(prhs[0]);
inVector2 = mxGetPr(prhs[1]);

/* get dimensions of the input matrix */
ncols = mxGetN(prhs[0]);

/* create the output matrix */
plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
/* get a pointer to the real data in the output matrix */
product = mxGetPr(plhs[0]);

/* call the computational routine */
dotProduct(inVector1, inVector2, product, (mwSize)ncols);
}



4- The main function that performs the main calculation, which is a dot product of two vectors in our example can be written as follows:

void dotProduct(double *M1, double *M2, double *M3, mwSize n)
{
mwSize i;
for (i=0;i<n;i++)
{
M3[0] += M1[i] * M2[i];
}
}


Here we have used mwSize instead of int which allows us to have cross-platform code.

5- Now, you can compile the code in MATLAB

mex multiplier.c


6- Now you can test the function in MATLAB. Remember that a MEX file can be directly executed in MATLAB:
[sourceode langaue=”matlab”]
A = [1 2 3];
B = [3 4 5];
multiplier(A,B)
ans =
26
[/sourcecode]

The whole sourcecode can be found as follows:

/*
* this is a very basic example of a matrix multiplier
*
*/

#include "mex.h"

void dotProduct(double *M1, double *M2, double *M3, mwSize n)
{
mwSize i;
for (i=0;i<n;i++)
{
M3[0] += M1[i] * M2[i];
}
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *inVector1;              /* input scalar */
double *inVector2;               /* 1xN input matrix */
size_t ncols;                   /* size of matrix */
double *product;              /* output matrix */

/* create a pointer to the real data in the input matrix  */
inVector1 = mxGetPr(prhs[0]);
inVector2 = mxGetPr(prhs[1]);

/* get dimensions of the input matrix */
ncols = mxGetN(prhs[0]);

/* create the output matrix */
plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
/* get a pointer to the real data in the output matrix */
product = mxGetPr(plhs[0]);

/* call the computational routine */
dotProduct(inVector1, inVector2, product, (mwSize)ncols);
}


## Investigating code performance in MATLAB – Part 1: For-loop

Here I want to show you a very simple example in MATLAB. Consider we have a large matrix (can be an image also) and we want to perform element-wise operation on the matrix elements. Let’s first create a large matrix:

A = rand(10000);


This creates a matrix with size: 10000 x 10000.

Consider we want to subtract each element of the matrix from the mean of the matrix. The mean of a matrix can be calculated using mean2 function. So the most straight forward technique is as follows:

tic
M = mean2(A);
for ii=1:10000
for jj=1:10000
A(ii,jj) = A(ii,jj) - M;
end
end
toc


Elapsed time is 12.007181 seconds.

Now, let’s try to increase the performance by removing one of the for-loops. The code would use a linear mapping between the element index and its row and column equivalent.

tic
M = mean2(A);
L = 10000^2;
for ii=1:L-1
kk = fix(ii/10000)+1;
jj = mod(ii-1,10000)+1;
A(kk,jj) = A(kk,jj) - M;
end
toc


Elapsed time is 19.745840 seconds.

Surprisingly, the execution time is increased. The reason is that in MATLAB, mod and fix are not operators but functions. And here by using 2 more functions we increase the execution time. The time is spent for calling those functions.

And the most efficient method works with removing both loops and using bsxfun from MATLAB instead.

tic
M = mean2(A);
C = bsxfun(@minus, A, M);
toc


Elapsed time is 0.721024 seconds.

So at the end we see that:

Two for-loops with no function call: 12.00 sec
One for-loop with two function call: 19:74 sec
No for-loop with bsxfun function: 0.72 sec

Posted in Linux, MATLAB, Optimization, programming, Software, Ubuntu | | 4 Comments

## Update status of your code output in one line (MATLAB)

Consider you want to report the status of your code in the command-line, but you do not want to print a new line each time. For instance, to show the percentage of the progress at the same line from 0%, to 10% and …. to 100%, you can do it in many possible ways, two of which I present here:

1- Using fprintf function to create dots for each step of the progress

fprintf('Processing data')
for ii=1:10
fprintf('.')
pause(0.5)
end
fprintf('\nDone!\n')


The result would be something like:

Processing data..........
Done!


2- Using fprintf function together with backspace operator to show the percentage of the progress:

fprintf('Processing data: \n');
for ii=10:10:100
fprintf('%d%%',ii)
pause(0.5);
fprintf('\b\b\b');
end
fprintf('\b\nDone.\n');


## Convert a grayscale image into a black-white image in MATLAB without im2bw

In MATLAB there is a function ‘im2bw’ which converts an input image into a black&white image. You can provide the function with a level of thresholding otherwise the default value is 0.5.

Here I show you how to implement the same functionality in a few lines of code. Check out the following function:

function S = my_im2bw(Ig,level)
S = Ig;
S(Ig > level) = 1;
S(Ig <= level) = 0;
end


Now let’s compare this function with the one of MATLAB:

Ig = rgb2gray(im2double(imread('test.png')));
level = 0.4;
B1 = im2bw(Ig,level);
B2 = my_im2bw(Ig,level);


The result is shown here: