Writing External C/C++ Functions for MATLAB

Anybody who has ever used MATLAB before will probably be aware that using loops is incredibly slow.  It is advised that you write as much of your code as possible using matrices and vectors.  However, if like me, your brain gets easily tangled into knots by multiple dimensions, or you want to do something particularly complicated, using matrices may not seem like a particularly attractive prospect.  I came across this problem (not for the first time) earlier today, and I decided I’d have a go at trying to solve it.

In these situations, execution time can be greatly reduced by writing and compiling C or FORTRAN functions, which you then call from MATLAB.  The form I will discuss here are C MEX files.  Simply by porting a loop from MATLAB to C, I managed to reduce the execution time of a simulation I was running from over 260 seconds, to about 20 (I’m working on getting it down to less than 10).

MATLAB Code

The original MATLAB code is below.  For those who are interested, it calculates the expected cost for buying a certain amount of energy at a given time, rather than waiting until closer to the time that energy is required for each discrete amount of energy one can buy.

for i = 1:length(amountToBuyInFM)
    cost(i) = amountToBuyInFM(i) * FMPrice
        + dot(consumptionDistribution(distroRange <= amountToBuyInFM(i)),
                    (distroRange(distroRange <= amountToBuyInFM(i)) - amountToBuyInFM(i)
            ).*(FMPrice-SBP_Delta))
        + dot(consumptionDistribution(distroRange >= amountToBuyInFM(i)),
                    (distroRange(distroRange >= amountToBuyInFM(i)) - amountToBuyInFM(i)).*(FMPrice+SSP_Delta)
            );
end;

This code is slow because MATLAB interprets the cost(i)= … function for every i in the range 1 to the length of amountToBuyInFM.  If we write this language in a compiled language it’ll be much faster, more fun, and probably a lot less ugly.

So… How do we code this in C++ in a way that MATLAB will understand?  Fairly easily, actually.  All you actually need is a function with the following signature.

#include "mex.h" // Don't worry about finding this file -- this is sorted out later.
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray * prhs[]);

Where, nlhs is the number of output parameters, plhs is a pointer to an array of the output parameters, nrhs is the number or input parameters, and prhs is a pointer to an array of input parameters.  This is called by MATLAB when your function is executed.  That’s pretty much it… There are just a couple of types, and functions that I’ve found quite useful, which I will discuss next.  By way of an example I’ll paste the code I wrote earlier (this is by no means the most efficient way of writing this code – but it is a LOT faster than running it in MATLAB…).  Those of a nervous disposition are welcome to skip over this, to the explanation that follows.

#include<stdio.h>
#include <gsl/gsl_cblas.h>
#include "mex.h"
#include "matrix.h"
#include <math.h>

#define DEBUG

#define AssertDouble( no ) if (! mxIsDouble (prhs[ no ] ) ) { mexErrMsgTxt( "Function 'MexGetCosts' requires parameter " #no " to be of type Double"); }

inline int find(const double * haystack, const double needle, const int haystacklength)
{
	for(int i =0; i < haystacklength; i++)
	{
		if(haystack[i] >= needle)
		{
			return i;
		}
	}
	return -1;
}

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray * prhs[])
{
	/* List of Arguments:
		1) Consumption Distribution (array)
		2) Distribution Range (array)
		3) Amount To Buy in Forward Market (array)
		4) Forward Market Price (scalar)
		5) Sell Delta (scalar)
		6) Buy Delta (scalar)
	*/

	//If the number of right hand side parameters is < 6 there is an error
	if(nrhs < 6)
	{
		mexErrMsgTxt("Function 'MexGetCosts' requires 6 arguments\n");
	}

	AssertDouble(0);
	AssertDouble(1);
	AssertDouble(2);
	AssertDouble(3);
	AssertDouble(4);
	AssertDouble(5);

	const double * consumptionDistro = mxGetPr( prhs[0]);
	const double * distroRange = mxGetPr( prhs[1]);
	const double * amountToBuyInFm = mxGetPr( prhs[2]);
	const double * fmPrice = mxGetPr( prhs[3]);
	const double * sDelta = mxGetPr( prhs[4]);
	const double * bDelta = mxGetPr( prhs[5]);

	const mwSize * distroRangeSize = mxGetDimensions(prhs[1]);
	const mwSize * fmAmountSize = mxGetDimensions(prhs[2]);

	//Creating the output matrix
	plhs[0] = mxCreateDoubleMatrix(1,fmAmountSize[1], mxREAL);
	//Give it a friendly alias
	double * output = mxGetPr(plhs[0]);

	//Now perform the calculation of the cost for each amount bought in the forward market
	int j;
	double acc = 0;
	double min = INFINITY;

	for(int i = 0; i < fmAmountSize[1]; i++)
	{
		//Get the first index in the distribution that gives us
		//the correct amount of power.
		j = find(distroRange, amountToBuyInFm[i],distroRangeSize[1]);

		if(j < 0)//Not found ...
		{
			output[i] = 0;
			continue;
		}

		output[i] = amountToBuyInFm[i] * fmPrice[0];
		acc = 0;

		//Dot product up to that amount (expected cost in wasted elec).
		for(int k = 0; k < j; k++)
		{
			acc += consumptionDistro[k] * (distroRange[k]-amountToBuyInFm[i]);
		}
		output[i] += acc* (fmPrice[0] - bDelta[0]);

		acc = 0;
		//Dot product over that amount (expected cost of late purchases).
		for(int k = j+1; k < distroRangeSize[1]; k++)
		{
			acc += consumptionDistro[k] * (distroRange[k] - amountToBuyInFm[i]);
		}
		output[i] += acc * (fmPrice[0] + sDelta[0]);
		if( output[i] < min)
		{
			min = output[i];
		}
	}

	//If we have more than one output argument, we'll output
	//the minimum in the second argument.
	if(nlhs > 1)
	{
		plhs[1] = mxCreateScalarDouble(min);
	}
}

Creating Matrices Etcetera

In order to create a double precision m-by-n matrix of real numbers, simply use the mxCreateDoubleMatrix(m, n, mxREAL), which returns an mxArray * pointing to your m-by-n matrix with all elements initialised to 0.  If you’d like to create a matrix of imaginary numbers, simply substitute mxREAL for mxCOMPLEX.  In order to obtain a pointer (double *) to the memory allocated to the real components of the matrix, you must use mxGetPr( mxArray * ), and for the imaginary components of complex matrices, mxGetPi.

In order to create a scalar value, simply use mxCreateScalarDouble, which takes a parameter stating what the value should be initialised to.  Easy.

Compiling

In UNIX this seems to be fairly easy too.  Simply use the ‘mex’ tool.  mex -help will give you a list of switches you can use, but I compiled the above code with the following command…

mex -lgslcblas mexcost.cpp

where mexcost.cpp is the filename of the source code in the previous box.

Nice and easy :D

Concluding Remarks

There’s a lot more to creating these functions.  For further information, see the Mathworks website.  A lot of mathematical functions have already been implemented by developers (surprisingly enough…).  In particular, check out BLAS, and Boost.